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ABSTRACT 


This  thesis  addresses  the  problem  of  computing  rapid  slew  maneuvers  for 
a  spacecraft  antenna  mounted  on  a  double-axis  gimbal  with  elastic  joints.  The 
performance  of  the  system  can  be  enhanced  by  designing  antenna  maneuvers  in 
which  the  flexible  effects  are  properly  constrained,  thus  reducing  the  load  on  the 
spacecraft  control  system.  The  motion  of  a  mass-spring-damper  system  is  shown 
to  be  analogous  to  a  spacecraft  antenna  slew  with  linear  dynamics.  This  model  is 
extended  to  a  nonlinear  double-gimbal  mechanism  with  flexible  joints,  which 
better  represents  real  spacecraft  antenna  dynamics.  Rather  than  increase 
maneuver  times  to  control  flexible  motion,  this  thesis  presents  optimal  solutions 
that  decrease  maneuver  times  while  allowing  designers  to  easily  constrain 
flexibility.  Since  it  is  impossible  to  recast  the  nonlinear  system  into  a  modal 
representation,  an  innovative  approach  is  used  to  map  the  nonlinear  dynamics 
into  a  linear  system  with  a  fictitious  force.  The  fictitious  force  captures  the  effects 
of  the  nonlinearities  so  the  vibrational  motion  can  be  constrained  for  a  time- 
optimal  slew.  It  is  shown  that  by  constructing  an  appropriate  optimal  control 
problem,  the  maneuver  time  for  a  flexible  DGM  can  be  decreased  by 
approximately  42%  compared  to  a  conventional  computed  torque  control 
solution. 
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I.  INTRODUCTION 


NASA  established  the  Tracking  and  Data  Relay  Satellite  (TDRS)  program 
in  1973  to  provide  continuous  global  services  to  critical  missions  in  low-earth 
orbit  while  minimizing  reliance  on  international  ground  stations.  Since  the  first 
launch  in  1983,  NASA  has  commissioned  and  operated  11  TDRS  in  space  (a 
twelfth  satellite  was  destroyed  in  the  Space  Shuttle  Challenger  explosion  in 
1986)  [1],  The  current  constellation  consists  of  nine  active  TDRS  at  various 
positions  in  geosynchronous  orbit  (see  Figure  1)  to  provide  global  relay  capability 
for  NASA  missions  such  as  the  International  Space  Station  and  the  Hubble 
Space  Telescope  [2], 
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Figure  1 .  TDRS  constellation,  from  [2], 


Each  TDRS  uses  multiple  antennas  to  service  its  customers  in  an  effort  to 
reduce  the  number  of  satellites  needed  on  orbit.  Third-generation  TDRS,  for 
example,  host  a  payload  with  two  single  access  antennas  and  a  multiple  access 
antenna  with  32  transmit  and  15  receive  elements  [2],  Each  of  these  antennas 


1 


must  be  steerable  in  order  to  service  moving  targets  in  the  low-earth  orbital 
regime.  This  presents  a  unique  problem  in  satellite  design  and  operation  in  that 
the  TDRS  attitude  control  system  must  continuously  accommodate  the  effects  of 
slewing  multiple  antennas  in  order  to  maintain  its  orientation  and  minimize  the 
associated  disturbances  on  the  spacecraft  bus  and  other  payloads.  The 
actuators  controlling  the  antenna  gimbals  must  supply  appropriate  levels  of 
torque  to  slew  the  antenna  at  a  desired  rate  in  addition  to  regulating  flexible 
effects  and  mitigating  pointing  errors  caused  by  the  motion  of  other  appendages. 
These  functions  are  accomplished  using  various  control  modes.  However,  for 
slews,  maneuvers  are  typically  done  slowly  to  avoid  exciting  the  flexible  modes  of 
the  satellite  structures.  In  contrast,  the  spacecraft  operators  endeavor  to 
maximize  the  utility  of  the  TDRS  by  servicing  as  many  customers  as  the  control 
system  will  allow.  Thus,  the  TDRS  program  faces  a  difficult  problem  in  that 
demand  for  access  is  high  and  availability  is  limited  due  to  the  capabilities  of  the 
various  satellite  control  systems  [3], 

A.  CONTROL  METHODS  FOR  A  FLEXIBLE  STRUCTURE 

Several  methods  exist  to  enhance  antenna  control  to  increase  user 
throughput.  The  antenna  can  be  modeled  as  a  manipulator  arm  connected  to  a 
controlling  body  through  joints  with  some  defined  range  of  motion.  The 
conventional  method  of  computed  torque  control  incorporates  a  simple  approach 
for  controlling  the  motion  of  such  an  arm  [4],  A  proportional-derivative  (PD)  or 
proportional-integral-derivative  (PID)  feedback  controller  is  common  in  computed 
torque  control  because  these  feedback  laws  can  provide  a  near-linear  response 
even  for  a  system  with  nonlinear  multi-body  dynamics.  A  simple  double  integrator 
can  then  be  used  to  model  the  motion  of  each  joint  of  the  arm  for  maneuver 
design.  The  feedback  mechanism  allows  the  controller  to  track  the  desired 
maneuver  profiles  [4],  [5],  This  method  of  simplified  maneuver  design  and 
implementation,  while  effective,  typically  dictates  that  the  system  be  operated 
conservatively  in  order  to  reduce  undesired  vibrations  in  the  antenna  and  other 
systems  on  a  satellite. 
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Another  method  for  antenna  control  applies  Pontryagin’s  principle  to 
develop  an  optimal  control  based  on  the  system  dynamics  and  given  boundary 
conditions  [6],  Pontryagin’s  principle  was  applied  in  [3]  and  [5]  to  find  minimum¬ 
time  solutions  to  the  TDRS  antenna-slewing  problem.  The  results  show  the 
potential  for  considerable  savings  in  terms  of  time  and  increased  customer 
access. 

For  simplicity,  most  models,  including  the  ones  in  [3],  assume  all  bodies 
are  rigid.  Realistically,  however,  space  structures  typically  exhibit  some  amount 
of  flexibility.  All  materials  possess  physical  properties  that  contribute  to  their 
response  to  forces,  specifically  impulses.  The  difficulty,  however,  lies  in  how  to 
model  this  flexibility.  Models  typically  require  comparison  to  actual  test  results  to 
ensure  accuracy.  As  pointed  out  in  [7],  mathematical  models  help  engineers 
understand  the  behavior  of  flexible  systems  in  response  to  (or  in  the  absence  of) 
external  input.  These  models,  however,  are  usually  simplified  to  avoid  nonlinear 
dynamics  because  the  nonlinearities  cause  unpredictable  behavior,  especially  in 
the  realm  of  flexible-body  motion.  Thus,  modeling  of  flexible  effects  in  mechanical 
systems  is  typically  linear.  This  thesis  explores  the  use  of  one  such  model  to 
optimize  antenna  maneuvers  while  taking  advantage  of  the  physical  properties  of 
materials. 

B.  THESIS  OBJECTIVES  AND  SCOPE 

A  significant  challenge  in  designing  control  systems  is  to  produce  rapid 
slew  maneuvers  while  minimizing  the  impact  of  flexible  motion.  An  antenna  must 
point  at  its  targets  with  a  certain  degree  of  accuracy  to  ensure  link  closure  and 
error-free  transmission.  A  rapid  slew  is  ineffective  if  it  causes  an  antenna  to 
“wobble”  off  target.  Figure  2  provides  a  simple  illustration  of  this  phenomenon. 
The  antenna  or  ground  station  must  increase  power  to  overcome  the  inaccuracy 
and  maintain  the  link.  To  develop  rapid  slew  maneuvers,  the  effects  of  flexibility 
must  be  considered.  To  this  end,  this  thesis  explores  the  concepts  of  modal 
analysis  by  looking  at  general  linear  mass-spring-damper  systems  with  two  and 
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three  degrees  of  freedom  (DOF).  These  linear  systems  are  notional 
approximations  of  two-axis  spacecraft  antennas  with  linear  dynamics.  Applying 
modal  analysis  to  these  linear  systems  helps  designers  to  understand  and 
constrain  vibrational  effects  of  flexible  systems. 


X 

target 


Figure  2.  Illustration  of  a  double-axis  gimbaled  antenna  vibrating  off 

target. 


Mapping  modal  analysis  concepts  to  the  antenna  slew  problem  is  difficult 
because  the  dynamics  of  the  antenna  motion  are  nonlinear.  Hence,  a  goal  of  this 
thesis  is  to  develop  an  approach  for  solving  rapid  antenna  slews  for  nonlinear 
antennas  that  enables  the  flexible  effects  to  be  properly  constrained  in  the  modal 
space. 

Chapter  II  describes  the  basic  concepts  in  modal  analysis  and  discusses 
the  value  of  frequency  and  impulse  response  functions  to  understanding  the 
behavior  of  the  flexible  system.  Two  methods  for  producing  a  modal  model  are 
presented,  one  in  the  frequency  domain  and  one  in  the  Laplace  domain. 
Chapter  III  expands  the  model  to  a  3-DOF  system,  representative  of  a  notional 
spacecraft  with  an  antenna  mounted  to  a  two-axis  gimbal.  Computed  torque 
control  is  used  to  implement  conventional  slews  to  illustrate  the  behavior  of 
flexible-body  motion  for  the  baseline  system.  In  Chapter  IV,  optimal  control 
theory  is  applied  to  minimize  slew  time  in  the  context  of  several  different 
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approaches  for  constraining  flexible  effects.  Finally,  in  Chapter  V,  these  results 
are  applied  to  a  more  realistic  spacecraft  antenna  model  with  nonlinear 
dynamics.  Chapter  V  presents  and  illustrates  an  innovative  approach  for 
marrying  modal  analysis  with  nonlinear  dynamical  systems  by  treating  the 
contributions  of  the  nonlinear  terms  as  a  fictitious  torque  input.  It  is  shown  that 
the  “new”  system  can  be  analyzed  and  optimized  as  a  series  of  linear  single-DOF 
systems  but  without  the  deleterious  effects  of  a  conventional  linearization 
process.  The  thesis  is  brought  to  a  close  with  some  summary  remarks  and 
suggestions  for  future  work  in  Chapter  VI. 
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II.  BACKGROUND  ON  MODAL  ANALYSIS 


In  modal  analysis  theory,  the  analysis  of  flexible  effects  is  conducted  in 
three  phases.  First,  one  must  derive  the  governing  equations  of  motion.  That  is, 
one  must  determine  the  values  that  populate  the  mass,  damping,  and  stiffness 
matrices  for  the  system.  These  second-order  models  are  applicable  because 
vibrational  amplitudes  are  usually  small  and  within  the  linear  range  of  material 
properties.  Second,  one  must  perform  a  free  vibration  analysis  of  the  system  by 
setting  any  external  forces  to  zero.  The  second  phase  of  analysis  produces  a  set 
of  natural  frequencies  and  damping  ratios  that  lead  to  a  matching  set  of  mode 
shape  vectors.  These  data  are  called  the  modal  properties  of  the  flexible 
system  [7],  The  mode  “shape”  is  defined  as  a  unique  “displacement  pattern”  [8] 
that  describes  the  “relative  displacements  of  all  parts  of  the  system”  [7],  Simply 
put,  a  vibration  “mode”  is  one  component  of  a  set  of  ways  by  which  a  system 
oscillates  in  response  to  an  input.  The  mode  shape  describes  how  each  part 
moves  relative  to  the  others.  Third,  the  forced  response  can  be  analyzed  to 
understand  the  behavior  of  a  flexible  structure  in  response  to  an  external 
input  [7],  which  is  of  interest  in  many  engineering  systems. 

The  general  form  of  the  governing  equation  of  motion  for  a  mass-spring- 
damper  system  with  n  degrees  of  freedom  is 

Mx(r) + Cx(r) + Kx(r)  =  f  (r)  (2.1) 

where  MgM"  is  a  diagonal  matrix  of  mass  values,  CeM'!X'!  is  a  matrix  of 
damping  coefficients,  KeK”  is  a  matrix  of  spring  constants,  f(?)el"!l  is  a 
vector  of  external  forces  acting  on  each  mass,  and  x(r)eR'!Xl  is  a  vector  of  the 
displacement  of  each  DOF  with  respect  to  its  equilibrium  position. 
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Consider  the  two-mass  system  shown  in  Figure  3. 


*i(Q 


m 


xx-0 


x2(t) 


m 


x2  =  0 


Figure  3.  Example  mass-spring-damper  system  with  two  degrees  of 

freedom. 


The  schematic  describes  a  2-DOF  system  in  motion  in  a  frictionless  environment. 
The  equations  of  motion  for  each  mass  are 


/j  =  mxxx  +  cx ,  +  kxx  -  cx2  -  kx2 


(2.2) 


f2  =  m2x2  -  cxx  -  kxx  +  cx2  +  kx2  .  (2.3) 

Using  (2.2)  and  (2.3),  the  mass,  damping,  and  stiffness  matrices  in  (2.1)  may  be 


assembled  as  M  = 


m , 

0 

r  -i 

i 

,  c  = 

c  — c 

0 

m2 

-c  c 

,  and  K  = 


k  -k 
-k  k 


The 


state  and  force  vectors  are  x(r)  = 


*i  (t) 
x2(t) 


and  f(0  = 


m 

m 


.  The  2-DOF 


system  is  used  in  this  chapter  to  illustrate  the  different  ways  in  which  the  modal 
properties  of  a  flexible  system  can  be  evaluated.  Undamped  motion  and  motion 
with  viscous  damping  will  be  studied  in  the  absence  and  presence  of  external 
forces  to  understand  the  behavior  of  the  system  shown  in  Figure  3. 


A.  NATURAL  MODAL  PROPERTIES 

A  system  without  external  excitation  will  exhibit  unique  modal  properties. 

These  properties  are  referred  to  as  natural  modal  properties  [7],  The  natural 

modal  properties  of  a  system  describe  the  behavior  of  a  forced  system  with 
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respect  to  certain  frequencies.  To  find  the  natural  modal  properties,  first  consider 
the  undamped  case  (c  =  0)  with  f(r)  =  0: 

Mx(0  +  Kx(f)  =  0.  (2.4) 


Assume  the  solution  to  (2.4)  is  of  the  form  x(t)  =  XeiM,  where  XeM"xl 
denotes  a  vector  of  time-independent  amplitudes.  This  leads  to  x(t)  = -co2Xeim , 
which  shows  that  the  system  oscillates  following  simple  harmonic  motion  with 
frequencies  co  and  amplitudes  X[7],  Substituting  x(t)  =  Xeim  into  (2.4)  gives  the 
general  eigensolution  of  the  form 

-M(02Xeim  +  KXei<0'  =  0 

(2.5) 

(K-ArM)X,=0 

where  Xr  is  the  r-th  eigenvalue.  The  corresponding  eigenvectors,  Xr ,  are 
scaled  versions  of  amplitudes  X .  For  a  system  with  n  degrees  of  freedom, 
r  —  and  X r=co2  where  cor  is  the  r-th  undamped  natural  frequency. 
Eigenvectors,  v,. ,  for  each  eigenvalue  must  be  computed  to  determine  the  mode 
shapes  of  the  system.  The  mode  shape  describes  how  the  system  responds  at 
the  corresponding  natural  frequency,  and  this  shape  is  dominant  at  that 
frequency  [9],  In  (2.5),  (K-ArM)=0  must  be  true;  otherwise,  X,  =0,  the  trivial 


solution.  Let  VeC'!Xl 

v  =  l 


be  a  matrix  whose  columns  are  eigenvectors  such  that 
I  y  ].  Thus,  the  modal  model  of  a  system  is 


completely  described  by  the  two  eigenmatrices 


and 


0 

0 

0 


0 

a2 

0 

0 


o  0 

o  0 

0 


(2.6) 


(2.7) 
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For  the  system  given  in  Figure  3,  it  is  desired  to  find  the  eigenvalues. 
Since  this  is  a  system  with  two  degrees  of  freedom  ( n  =  2)  there  are  two  modes, 
two  eigenvalues,  and  two  eigenvectors.  Applying  (2.5)  gives 


k  -  Arm, 
-k 


k-X 


/A, 


X  =0 


(2.8) 


The  eigenvalues,  A,. ,  may  be  solved  as 


det 


k-X, 

-k 


m. 


k  -  Xm2 


=  mxm2Xr  -  (ml  +  m2  )kXr  =  0 


(2.9) 


Aj  =(o{2  =0 
1  „2  (ml  +  m2)k 

A2  CiJ  2 


(2.10) 


The  eigenvectors  can  be  found  by  substituting  the  eigenvalues  back 
into  (2.8)  and  solving  for  \r  (MATLAB’s  eig  command  can  do  this  symbolically). 

For  A,  =  0, 


k-\mx  -k 
-k  k-X1m2 


1 

1 

3 

1 

0 

k  -k 

0 

r  *  -i 

0 

-k  k  —  X]m2 

0 

-k  k 

0 

o 

o 

0 

v„-v21=0  (2.11) 


vn 

vu 

a 

V2 1 

=  a 

vu 

where  a  is  a  scaling  factor.  Choosing  a= —  gives 
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(2.12) 


For  A2  = 


(mt  +m2)k 
mxm2 
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Choosing  a- — , 


22 


Therefore, 


v2 


ra, 


1  1 


(2.15) 


(2.16) 


In  the  first  mode  (described  by  the  first  column  of  V),  both  masses  move  in 
the  same  direction  with  the  same  amplitude,  since  both  elements  of  v,  have  the 
same  sign  and  value.  The  second  mode  is  characterized  by  the  masses  moving 
in  opposite  directions  due  to  the  negative  sign,  and  the  relative  amplitude  of 

YYl 

motion  is  determined  by  the  ratio  of  the  masses  (i.e.,  —  [10]).  Because  the 

m] 

columns  of  (2.16)  are  eigenvectors,  they  are  scalable  in  amplitude.  However,  the 
shapes  of  the  modes,  and  hence  the  relative  behavior  of  the  two  masses,  will  not 
change  [7],  The  undamped  free  motion  of  the  2-DOF  system  is  therefore 
completely  described  by 


[K]  =  [c°r]  = 


0  0 
0  (ni\+m2)k 

m1m2 


and 


V  = 


1 


1 


m , 
1 


(2.17) 


where  [Ar]  denotes  a  diagonal  matrix  of  eigenvalues  and  [ar]  denotes  a 
diagonal  matrix  of  natural  frequencies,  and  r  =  l...n. 
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B.  DIAGONALIZING  THE  MASS  AND  STIFFNESS  MATRICES 


An  important  property  called  orthogonality  can  be  used  to  diagonalize  the 
mass  and  stiffness  matrices  in  order  to  transform  the  original  spatial  model, 
which  describes  the  motion  in  terms  of  displacements,  x,  to  the  so-called  modal 
model,  which  describes  the  motion  in  a  new  set  of  coordinates.  The 
transformation  relies  on  the  fact  that 


M  =  VrMV 
K  =  VrKV 


(2.18) 


where  Mel"  is  a  diagonal  matrix  of  modal  masses,  mr ,  and  KgI“  is  a 

_ 

diagonal  matrix  of  modal  stiffness  values,  kr ,  leading  to  Ar  =  co;  =  —  .  The  modal 

mr 

mass  and  modal  stiffness  matrices  are  not  unique.  Rather,  they  are  based  on 
how  the  eigenvectors  are  scaled.  However,  the  ratio  of  these  quantities  is  unique 
and  fixed,  giving  unique  modal  frequencies,  as  in  (2.17).  The  modal  mass  matrix 
is  generally  used  to  normalize  the  mode  shape  vectors.  The  normalized  mode 
shape  vectors  can  then  be  used  to  derive  useful  quantities  such  as  effective 
mass  or  effective  stiffness  [7],  The  mass-normalized  mode  shape  vector  (pr  eR'lxl 
is  calculated  as 


0,-  = 


(2.19) 


such  that  column  matrix  OgR™  is 


o  = 


0i 


02 


(2.20) 


which  has  the  property 

dr'^O7)-1  (2.21) 

For  the  2-DOF  system  given  in  Figure  3,  the  modal  masses  and  stiffness 
values  are 
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and 


1 

i 

r 

1 

,  m2 

_ 

m. 

0 

1 - 1 

M  =  yrMV  = 

m2 

i 

m. 

— 

0 

m2 

1 

mx 

-1 

1  1 

1 

1 

m2 

m. 

w2 

ml 

—m2 

m2 

1 

1 

m[  +  m2  0 


m{ 


(2.22) 


1  1 

,  m2 

K  = 

VrKV  = 

m2 

k 

k 

1 - 1 

m. 

- 1  1 

-k 

c 

l 

m. 

1  1 

0 

0  1 

1 

2_ 

= 

1(ml  +  m2)  J(ml  +  m2) 

/C  rv 

i 

mx 

mx 

mx  , 

i 

1 

0  0 

ml 


(2.23) 
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Using  (2.22)  and  (2.23),  it  is  possible  to  verify  that  [Ar]  =  M  'k  as 
required: 


mx  +  m0 

0 

-1 

0 

0 

[A,.]  =  M  K  = 

0 

2 

ffl, 

— -  +  m2 

0 

k{mx  +  m2 )” 

mx  +  m2 


m , 


m2(m,  +  m2) 


0  0 


k[mx  +m2) 


m. 


0  0 
0  (ni\+m2)k 

m1m2 


(2.24) 


The  mass-normalized  mode  shape  vectors  are 


o  = 


01  02 


(m.  +  m2 ) 


■i/  m. 


m. 


m. 


N  1/ 
\/2 


^  (m,  +  m2  )m2  J 


(m,  +  m2)  /2 


in. 


X  1/ 

\/2 


(/».,  +  m2)m2  J 


(2.25) 


Figure  4  and  Figure  5  show  an  example  of  the  relationship  between  the 
mode  shapes  and  the  displacement  of  each  mass,  for  a  sample  system  where 


m1  =  m2  =  k  =  l.  The  values  of 


V2  V2 

2  2 


0.7071  0.7071  J  (solid 


line)  mean  that  each  mass  will  move  in  phase  and  with  the  same  amplitude  at  the 
frequency  of  the  first  mode,  as  seen  in  Figure  5(a).  For  the  second  mode,  each 
mass  will  move  in  opposite  directions  but  with  the  same  magnitude.  The  motion 
is  out  of  phase.  In  Figure  5(b),  which  shows  the  free  motion  of  the  2-DOF 
system,  mass  1  oscillates  at  the  same  amplitude  as  mass  2  but  in  the  opposite 
direction. 
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Amplitude 


2 

1.5 

1 

0.5 

0 

-0.5 

-1 

-1.5 

-2 

Figure  4.  Mass-normalized  mode  shapes  for  a  2-DOF  system. 


Mass  # 


0  2  4  6  8  10 

(a) 

Time  (sec) 


Time  (sec) 


Figure  5.  Normalized  displacements:  (a)  mode  1;  (b)  mode  2. 
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C.  UNDAMPED  CASE  WITH  AN  EXTERNAL  FORCE 

Adding  a  forcing  function,  f  (t) ,  to  the  2-DOF  system  gives 

Mx(0+Kx(0=f(0-  (2.26) 

The  presence  of  the  forcing  function  changes  the  procedure  for  modal  analysis. 

In  this  case,  the  force  vector  can  be  written  as  f  =  Fem  where  F  e  C"xl  is  a  vector 
of  complex  amplitudes  and  co  is  the  frequency  [7],  The  variable  co  does  not 
represent  an  output  frequency  as  in  (2.5);  rather,  it  denotes  an  input  frequency 
for  the  forcing  function.  The  eigensolution  in  (2.5)  then  becomes 

(K  - co2M)Xeim  -  Fel0>t  .  (2.27) 

Following  the  procedure  in  [7],  let 

X  =  H(to)F  (2.28) 


such  that  H  e  M',x"  and 

H(©)  =  (K-6)2M)-1.  (2.29) 

In  (2.29),  matrix  H  is  called  a  receptance  frequency  response  function 
(FRF);  it  represents  the  system’s  displacement  in  response  to  an  input  given  in 
the  frequency  domain  [7],  Using  the  modal  properties  of  the  system,  we  may 
multiply  by  the  matrix  <t>  and  its  transpose,  as  shown: 

®rH(«)  =  Or(K-arM)®  =  OrK®-®WMO  .  (2.30) 

Since  <t>  is  mass-normalized,  the  quantity  OrKO  gives  a  diagonal  matrix  of 

k 

-zr  =  (o;  values.  It  also  follows  that  =  «2®7M®  =  diag(ro2). 

mr 

Therefore,  (2.30)  may  be  rewritten  as 

OrH(«)  10  =  diag[(w2-«2)]  (2.31) 

and  further  manipulated  to  give 

H(ffl)"1  =(0rrIdiag[(©2-6)2)]0-1  .  (2.32) 
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Recalling  that  (AB)  1  =  B  'A  1 ,  we  obtain 


H(«)  =  0(diag[(«;-(y2)]) .  (2.33) 

From  (2.33),  it  is  apparent  that  matrix  H  is  symmetric.  Because  H 
represents  the  relationship  between  displacement  and  force,  its  components  can 
be  computed  as 

H*(®)  =  fr  (2-34) 

F, 


where  j  =  \...n  and  k  =  \...n.  In  other  words,  Hjk  is  the  FRF  for  the 
displacement  of  the  j'-th  mass  due  to  the  k-\h  element  of  the  force  vector, 
superposed  over  all  the  modes.  The  same  relationship  will  be  explored  later  in 
the  Laplace  domain.  This  “element-wise”  relationship  is  useful  to  study  the 
response  due  to  certain  forces  on  specific  masses  in  the  system  [7],  Equation 
(2.33)  could  also  be  used  to  find  the  full  matrix  of  FRFs. 


H  ;,(©) 


I 


<E>  O, 

jr  kr 


co: 


CO 


"  V  V  "  R 

y  Y  jr  ykr  _  y  r^jk 

hm,\co;-co2)  “l  co2  -co2 


(2.35) 


The  modal  constant  rR.teClxl  is  the  specific  receptance  linking  coordinates  j 

and  k  for  mode  r  [7],  The  notion  of  the  modal  constant,  or  “residue,”  will  also  be 
seen  later  in  the  Laplace  domain. 

Equation  (2.35)  highlights  the  principle  of  superposition  with  respect  to 
modes — the  total  response  of  a  system  is  the  sum  of  the  responses  of  each 
mode.  In  other  words, 


H(g»  =  X,H«o) 

r= 1 


where  ,  H(o>)  denotes  the  FRF  for  mode  r .  In  other  words, 


(2.36) 


n  n 

*  II » ►  • 


r= 1  k=\ 


(2.37) 


18 


For  the  2-DOF  system  shown  in  Figure  3,  the  first  mode  dictates  motion  of 
both  masses  in  the  same  direction  and  with  the  same  amplitude.  Thus,  the 
receptance  FRFs  are  the  same. 


(m]  +m2 )  /2(ml  +  m2) 
-co2 


-1 


or  (m,  +  m2) 


Using  (2.38)  gives 


|Hn 

h21 


!H12 

IT 

l'n22 


-1  -1 

co2(ml  +  m2)  0J2(mt  +  m2) 
-1  -1 

co2(ml  +  m2)  or(m]  +  m2) 


(2.38) 


(2.39) 


For  the  second  mode,  the  amplitude  of  motion  is  the  same  for  each  mass, 
but  the  directions  oppose  each  other.  The  symbolic  representation  for  the  second 
mode  is  cumbersome,  so  0H(ro)  is  not  provided  here.  A  numerical  example  is 

given  later  to  clarify  the  details. 


D.  VISCOUS  DAMPING 

So  far,  only  the  undamped  case  has  been  studied.  More  realistic  is  the 
presence  of  damping  in  some  form.  Two  general  cases  are  discussed  in  [7]: 
structural  (or  hysteretic)  damping  and  viscous  damping.  While  realistic  to  some 
degree,  there  are  advantages  and  disadvantages  for  each  case.  In  this  thesis, 
only  viscous  damping  is  considered. 

In  the  case  of  viscous  damping,  we  first  seek  to  determine  the  natural 
modes  of  the  system  by  examining  the  unforced  system,  for  which  the  general 
equation  of  motion  is 

Mx(f )  +  Cx(f )  +  Kx(f )  =  0 .  (2.40) 

The  properties  of  each  matrix  are  the  same  as  previously  defined  for  an  n  -DOF 
system.  For  the  unforced  case,  with  a  solution  of  the  form  x  =  \eh ,  where  A  =  ico , 
the  equation  of  motion  yields 
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(MA2  +  CA  +  K)X  =  0  . 


(2.41) 


The  solution  to  this  equation  is  a  complex  eigenproblem,  with  complex- 
conjugate-pairs  of  eigenvalues  A  and  eigenvectors  v  such  that 

(MA2+CA  +  K)v  =  0  .  (2.42) 

The  bar  symbol  over  the  eigenvalue  and  eigenvector  variables  is  used  to 
distinguish  them  from  the  ones  computed  for  the  undamped  free  case.  Each 

eigenvalue,  A,. ,  can  be  expressed  as 

A,  =  +  i-vA-fT)  (2.43) 

so  that  the  eigensolution  of  (2.42)  can  be  described  as 

K’K 

and  ,  (2.44) 

vr,vr* 

where  r  =  l...n  .  In  (2.44),  the  *  notation  is  used  to  denote  the  complex  conjugate, 
(i.e.,  X*r  =  cor - iyjl - £r2 ) ).  The  orthogonality  conditions  are  much  more 
complicated  than  in  the  previous  case.  Let  v  be  the  eigenvector  for  the 
q-\h  mode;  the  same  notation  applies  for  A  .  Ewins  [7]  derives  the  following 
orthogonality  conditions: 

(Ar  +  A9)yjMyr  +  v9Cvr  =  0 

and  (2.45) 

A;.AXMv,.  +  v9Kv,.  =  0 
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If,  and  only  if,  v?  =  v*  and  Xq  =  ft),.  |-£r  -iyj  1  -  £,2  j  =  A*  then  it  is  possible  to 

determine  the  so-called  natural  frequencies  cor  and  critical  damping  ratios  of 
the  system.1  In  this  case,  (2.45)  leads  to 


v  "Cv 

2a)£r  = 

v  Mv. 


m.. 


CO, 


and 

Kv, 

m. 


v  Hxr- 

^7 


Mv. 


(2.46) 


where  the  superscript  H  denotes  the  Hermitian  (complex  conjugate)  transpose. 
The  similarity  to  the  conditions  in  (2.18)  is  evident. 

The  presence  of  damping  makes  the  analysis  of  the  forced  system  given 
in  (2.47)  more  difficult  to  manipulate  than  for  the  undamped  case. 

Mx(r)+Cx(r)+Kx(r)  =  f(r)  .  (2.47) 

To  develop  the  FRF  for  the  forced  response,  first  recast  (2.47)  into  state-space 
form  [1 1],  This  can  be  accomplished  using  the  following  procedure  if,  and  only  if, 
the  roots  of  the  characteristic  equation  of  the  system  have  multiplicity  one.  Define 
state  variable  yeM2"xl  as 


y  = 


X 

X 


(2.48) 


so  that  the  unforced  equation  of  motion  can  be  recast  as 

I  C  M  ]y  +  [  K  0  ]y  =  0  .  (2.49) 

This  produces  a  system  of  n  equations  with  2 n  unknowns,  which  is  not  helpful. 
Therefore,  we  make  use  of  the  identity  equation 

I  M  0  ]y  +  [  0  -M  ]y  =  0  (2.50) 

1  This  “natural  frequency”  matches  the  “undamped  natural  frequency”  for  single-DOF 
systems  or  systems  exhibiting  proportional  viscous  damping  only  [7]. 
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together  with  (2.49)  to  form  a  set  of  2 n  equations  to  give 


c 

M 

y  + 

K 

0 

M 

0 

0 

-M 

This  equation  simplifies  to  canonical  state-space  form 

Ay  +  By  =  0 


(2.51) 


(2.52) 


where  AgR2”2"  and  BeR2'1*2'1  are  each  square  matrices.  The  eigensolution, 
assuming  a  solution  to  (2.52)  of  the  form  y  =  Y<?Ir,  produces  2 n  complex- 
conjugate-pair  eigenvalues  and  2  n  complex-conjugate-pair  eigenvectors 
(pr  eC2”xl .  The  eigenvalues  and  eigenvectors  satisfy  the  general  eigensolution 

(X,.A  +  B)(pr  =  0  (2.53) 


where  r  =  l...2n.  Orthogonality  is  useful  here  as  well.  Let  i i)  e  C2nx2n  be  a  matrix 
whose  columns  are  the  eigenvectors,  (pr ,  of  the  uncoupled  system.  Then,  we 
may  write 


\ar]  =  $TA& 
[br\  = 


(2.54) 


such  that 


(2.55) 


where  ar  and  br  are  the  r-th  elements  of  diagonal  matrices  [a,.] eC2”12" 
and  [z?,.]  e  C2nx2" ,  respectively.  It  is  important  to  note  that  the  transpose  operation 

in  (2.54)  is  the  non-conjugate  transpose  of  the  eigenvector  matrix  (versus  the 
Hermitian  transpose  used  in  (2.46))  [7], 

The  forcing  function  f  can  now  be  expressed  as 


P  = 


F 

0 


(2.56) 
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where  PeC2"xl.  To  find  the  FRF,  recall  the  procedure  used  in  (2.28)-(2.35).  The 
FRF  can  then  be  written  as 


X 

icoX 


_y  <P^<Pr_ 
r= i  Ur(i(£  ~  A.r) 


(2.57) 


which  leads  to  the  following  individual  receptance  FRFs  for  the  forced  system 
(with  viscous  damping): 


H*(0» 


s 


f 

$j,Ar 

,  KK 

A 

^ar[co£r+i(o)-CQrJl-C; 

'))  «r(®rCr+'(®  +  ®r>/1-C  i 

0), 

.(2.58) 


The  relationship  in  (2.55)  points  to  a  useful  coordinate  transformation 
provided  by  the  eigenvector  matrix  i3 .  Let  y  =  $q  where  q  e  C2"xl  is  the  so- 
called  modal  coordinate.  From  (2.53),  one  can  see  that  q  =  /'.  Asa  result,  (2.52) 
becomes 


A#q  +  B#q  =  P 


(2.59) 


The  transformation  by  d  puts  the  equations  of  motion  in  terms  of  modal 
coordinates  to  simplify  further  analysis.  Premultiplying  by  A  yields 

tfrAtfq  +  tfrBtfq=tfrP  .  (2.60) 


So,  for  the  r-th  mode,  (2.54)  can  be  used  to  get 

tfrrAtf,.q,.  +  tf,rBtf,.q,.  =  ifrV  =  P,' 
arq  +  ^q,.  =  P; 


(2.61) 


Here,  the  r-th  mode  corresponds  to  the  r-th  element  in  the  vector  q  [11],  This 
transformation  will  be  discussed  further  below. 
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E.  PROPORTIONAL  VISCOUS  DAMPING 

A  special  case  of  viscous  damping  is  called  proportional  (or  classical) 
damping  [12].  Proportional  damping  offers  advantages  in  analysis  over  the  more 
general  case  in  that  the  mode  shapes  of  a  proportionally  damped  system  are  the 
same  as  those  for  the  undamped  system,  and  the  natural  frequencies  are 
similar  [7],  The  general  equation  of  motion  is  the  same  as  in  (2.40),  but  now  the 
damping  matrix  is  proportional  to  the  mass  and  stiffness  matrices  such  that 

C  =  /5K  +  yM  (2.62) 

where  /3  and  y  are  scalar  constants.  As  a  result,  C  is  diagonalizable  by  the 
eigenvector  matrix  V  to  produce 

C  =  VrCV  (2.63) 

where  the  r-th  diagonal  element  of  CgR'™  is  the  modal  damping  constant,  cr, 
for  the  r-th  mode. 

For  the  2-DOF  system  in  Figure  3,  the  modal  damping  matrix  is 


1 

1 

1 

C 

-c 

-c 

c 

1 

777, 

777, 

777, 

1 

1 

0  0 

=  (m1  +  m2)  (, mx+m2 ) 

in  |  ni\ 

0  0 
=  o  c(”h*m1f 

777, 

k 

because  (5  =  -  exists  and  7  =  0. 

c 

The  modal  transformation  of  a  proportionally  damped  system  is  also  more 
straightforward.  The  unforced  system,  from  (2.18)  and  (2.63),  is 


l  _"T 

777, 

1  1 


(2.64) 
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(2.65) 


Mq  +  Cq  +  Kq  =  0 

where  q  =  V _1x .  Thus,  the  forced  case  is  written  in  modal  coordinates  as 


Mq  +  Cq  +  Kq  =  Vrf (f )  =  f  '(f)  ■ 


(2.66) 


For  the  r-th  mode, 


mAr  +  crqr  +  krqr=fr\t) 


(2.67) 


Equation  (2.67)  is  clearly  the  equation  of  motion  for  a  set  of  single-DOF  mass¬ 
spring-damper  systems,  which  represents  the  modes  of  a  multiple-DOF  system. 
Assuming  a  solution  to  the  equation  of  the  form  q  =  QeAf ,  where  A  is  complex 
due  to  the  presence  of  the  damping  terms,  the  r-th  complex  conjugate 
eigenvalue  is  computed  from  the  characteristic  polynomial 


mrX2r  +  cr\  +  kr=  0 


A  =  ■ 


■  ±4 


—2 

C. 


■  4m  k,. 


2  m,. 


(2.68) 


where,  because  of  the  proportional  relationship  in  (2.62), 


and 


P<Dr  ,  Y 

2  2  cor 


The  damped  natural  frequency  co'  is  written  as 


or  =  a. 


We  note  that  as  long  as  is  small,  co'  ~  cor . 


(2.69) 


(2.70) 


For  the  2-DOF  system  shown  in  Figure  3,  by  inspection, 
co1  =  co[  =  =  Aj  =  0 ,  and  for  the  second  mode, 


25 


k(mx  +  m2 ) 


(2.71) 


co2  = 


1  k2 

k(ml  +  m2  )7 

2 

ml 

2 

m7 

— -  +  m2 

mx 

C2  = 


c(m,  +  m2y 

2 

in. 


c(m ,  +  m2 ) 


2  yfnA 


If  2 

\ 

f 

m7 

+  m2 

'  V  wi 

J 

V 

k(m]  +  m2  Y 

2 

m, 


\ 

2w  j  1 

C 

r  2  \ 

— -  + 

J  V 

l  Wl  J 

(2.72) 


For  analysis  of  the  forced  response,  it  is  possible  to  make  use  of  (2.29)  to 


get 


H(ry)  =  (K  +  icoC  -  arM)  ' 

VrH((y)V  =  Vr  (K  +  icoC  -  (y2M)  '  V  . 

Thus,  the  individual  FRFs  are  represented  by 


V  V 

jr  kr 


kr  +  icocr 


( o2mr 


(2.73) 

(2.74) 


(2.75) 


where  each  term  in  the  sum  is  the  FRF  corresponding  to  the  r-th  mode. 

For  the  2-DOF  system  shown  in  Figure  3,  the  same  symmetry  exists  in  H 
when  assuming  proportional  damping.  For  the  first  mode,  r  =  1,  the  FRFs  are  the 
same  as  for  the  undamped  system  (this  may  not  always  be  the  case)  in  (2.39). 

For  r  =  2 , 


2H„(ro) 


V  V 

M2  y  12 


k2  -  iCOc2  -  co  m2 


(  V 

V  mx  J 


k  ( m,  +  m, )  c(m.+  m1 )  , 

V  +iCO  V  ~(0~ 


m: 


m. 


m: 


f  ivy 
V  mi 


+  m, 


k(m[  +  m2 )  +  icOc ( m ,  +  m2 )  -  ft)  ( m, m2  +  mxm2 ) 
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(2.76) 


and  so  on.  Numerical  examples  are  discussed  below  that  further  illustrate  the 
features  of  the  modal  analysis. 


F.  NUMERICAL  EXAMPLE 

Assume  a  nominal  case  where  ml=m2=c  =  k  =  1 .  For  the  undamped  case 
with  free  vibration, 


w= 


0 

0 


0 

( ml  +  m2 )  k 
mlm2 


0  0 
0  2 


(2.77) 


,  m2 

1 - - 

1  -1 

m. 

= 

1 

1  1 

1  1 

It  follows  that  ft),  =  0  and  co2  =  J 2.  Column  vector  v2  indicates  that  the 
amplitudes  of  motion  of  both  masses  in  mode  2  will  be  equal  but  in  opposite 
directions.  Thus,  the  modal  matrices  and  mass-normalized  mode  shape  vectors 
are 


with 


ml  +  m2 
0 


0 


+  m2 


2  0 
0  2 


k(ni\  +  m2 ) 


m. 


0  0 
0  4 


(2.78) 


27 


(2.79) 


o  = 


(»?,  +m2)~ 

(m,  +m2) 


m, 


N  1/ 
\/2 


( ml  +  m2)m 


m. 


2  J 

\V2 


(m,  +  m2 )  m 


2  y 


V2  V2 
2  2 
V2  V2 
2  2 


Note  that  with  different  scaling  of  the  eigenvectors  in  V ,  the  modal  matrices  will 

contain  different  values,  but  the  relationship  Ar  =  co2  =  —  will  be  preserved. 

mr 

To  analyze  the  frequency  response  in  an  undamped  system  with  an 
external  force,  use  (2.33)  to  get 


H(ffl)  = 


V2  _ 

_V2 

-i 

V2 

V2 

2 

2 

-co2 

0 

2 

2 

V2 

V2 

0 

2 -(o2 

2 

2 

2 

2 

1 

1 

1 

2  co2 

2  (or 

!-2) 

or  (or 

- 

2) 

1 

1 

1 

(O2 

(or -2) 

2  co2  2( 

CO 

2 -2) 

(2.80) 


Alternatively,  each  element  of  H  can  be  computed  individually  using  (2.35).  For 
example, 


Hi2(g» 


'vr 

'V2' 

'  ■&' 

'vr 

2  J 

2  J 

2  J 

2  J 

2 -ry2 


1 

co2(co2  -2) 


(2.81) 


A  plot  of  the  impulse  response  for  the  FRF  in  (2.81)  is  shown  in 
Figure  6(a)  as  an  oscillating  linear  growth  toward  infinity.  Response  to  other 
inputs  would  show  growth  to  infinity  as  well.  This  is  not  helpful  in  studying  the 
oscillatory  effects  only,  because  the  rigid-body  motion  of  the  system  and  the 
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oscillations  are  superposed,  (i.e. ,  the  time  history  of  Figure  6(a)  is  the  sum  of  the 
rigid  motion  in  Figure  6(b)  and  the  oscillatory  motion  in  Figure  6(c)). 


Time  (sec) 

Figure  6.  Undamped  impulse  responses  in  time  domain: 
(a)  H12  —  jH12  +  2H12 ,  (b)  j Hj2  ,  (c)  2H12  ■ 


For  reference,  the  FRF  for  the  first  mode  (from  (2.81 ))  is 

I 

,H12(£0)  =  -4  (2.82) 

co~ 

and  the  FRF  for  the  second  mode  is 
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1 

2 


(2.83) 


2H12(ft>)  - 

2-oj 


Since  the  damping  matrix  of  the  system  shown  in  Figure  3  is  proportional 
to  the  stiffness  and  mass  matrices  (/J  =  l,  7  =  0),  this  system  exhibits 
proportional  damping  and  is  straightforward  to  analyze.  However,  for 
completeness,  the  general  complex  modal  case  is  computed  to  illustrate  the 
steps.  First,  the  unforced  case  is  analyzed  to  compute  the  natural  frequencies 
and  damping  ratios.  First,  we  find  the  eigenvalues  in  (2.42). 


det(l2M  + AC  +  k)  =  0 


f 

- 

- 

- 

- 

r  -A 

A2 

0 

A 

-A 

1  -1 

+ 

+ 

V 

0 

A2 

-A 

A 

-1  1 

L 

j  A~+A  +  l  —  A  —  1 
=  det  _  _  _ 

— A  - 1  A’  +  A+l 

=  Ar4  +  2Ar3  +  2Ar2  =  0 
A,4  +  2A,3  +  2  A,.2  =  Ar2  (A,.2  +  2Ar  +  2)  =  0 


(2.84) 


Therefore,  Xi=X*=0.  The  quadratic  formula  yields  A2=-l±/  from  the  second 
term.  Following  the  same  process  in  (2.11),  for  Xl=X1*=0,  we  obtain  the  first 
eigenvector: 


A,  +  Aj  + 1 

-A,  - 1 

0 

1 

-1 

0 

1  -1 

0 

-A,  -1 

A|”  +  Aj  + 1 

0 

-1 

1 

0 

o 

o 

0 

(2.85) 


For  A2  =  —1  +  i , 
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(—1  +  /)  +  (— l  +  z')  +  l 

-(-t  +  z)-1 

0 

-(-t  +  z)-1 

(— 1  +  z)  +  (— l  +  z')  +  l 

0 

— z  -i 

0 

1 

1 

0 

-i  -i 

0 

0 

0 

0 

v2  = 

1 

±i 

0 

z 

-1 

0 

(2.86) 


Recall  that  the  eigenvectors  occur  in  complex  conjugate  pairs,  hence  the  ±  sign. 

Use  these  eigenvectors  and  (2.46)  to  find  the  natural  frequencies  and 
critical  damping  ratios  for  each  mode.  For  mode  1 , 


c o , 
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Vj^MVj 
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t 

0 

t 

0 

1 

l 
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=  0 


(2.87) 


which  follows  since  A,  =  V  =0 .  For  mode  2, 
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_  v2hCv2  _  ^ 

1 

->] 

1 

-1 

-1 

1 

1 

-1 

ni2 

v2"Mv2 

[> 

-1 

1 

[  0 

0 

1 

1 

-1 

_ 

_  v2HKv2  _  ^ 

1 

-•]’ 

1 

-1 

-1 

1 

r  i 

-i 

m2 
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=  1  =  2 
2 


=  2 


(2.88) 


Therefore,  co2=J 2  and  £2  =  .  As  a  check,  we  may  verify  that  these 

values  reproduce  the  eigenvalues  given  by  (2.84): 
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i,=w2(-t;2+i^Q)=sl2 


V 


■Ji  .1,  fVP 

- +ljl-  - 

V  2  J 


-\-\-i 


(2.89) 


X2  =  -1  —  i 


Now  we  wish  to  derive  the  uncoupled  modal  motion  of  the  system.  The 
state-space  form  of  this  system  (from  (2.51 ))  is 


1 

-1 

1 

0 

1 

-1 

0 

0 

-1 

1 

0 

1 

y+ 

-1 

1 

0 

0 

1 

0 

0 

0 

0 

0 

-1 

0 

0 

1 

0 

0 

0 

0 

0 

-1 

(2.90) 


Using  the  eigensolution  form  given  in  (2.53),  the  eigenvalues  may  be 
computed.  The  same  eigenvalues  as  computed  above  should  be  obtained.  The 
results  are 


det[l,.A  +  B]  =  0 
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=  V  +  2V  +  2V 
=  Xr2  (Xr2  +  2X,.  +  2)  =  0 


(2.91) 


Again,  ^  =0  and  A2  =-l±i .  The  eigenvectors,  however,  will  be  different 
because  the  system  is  now  in  the  uncoupled  state-space  form.  For  X,  =  0 , 
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(2.93) 


Recall  that  the  eigenvectors  occur  in  complex-conjugate  pairs,  which  is  the  case 


obtained  in  (2.93).  The  full  matrix  d  is  of  the  form 


(Pi  Vl  <Pl  (Pi 


and  produces  matrix  [ar]eC2nx2n  of  the  form  diag 


*  * 


Cly  CI2  Cly  0-2 


and 


\b,.\  =  diag 


bx  b2  b*  b2 


.  From  (2.54), 
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[ar]  =  tfrAtf  = 


r 

1 

1 

1 

-1 

1 

0 

1 

1  1 . 

0  0 

-1  1 

1  _  1  . 

-1 

1 

0 

1 

—  ±—l 

-\-  l 

1 

0 

0 

0 

2  2 

2  2 

0 

1 

0 

0 

-1 

0  0 

0  -2 


i  i±i/ 

2  2 

i  -ITI, 

2  2 


0 

0 


-1 

1 


(2.94) 


r 
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1 

-1 

0 

0 

1 

1  1 

0  0 

-1  1 

1  1 

-1 

1 

0 

0 

—  ±—  I 

—  +—i 

0 

0 

-1 

0 

2  2 

2  2 

0 

0 

0 

-1 

0 

0 


0 

-2  +  2  i 


1 

1 

0 

0 


1  ±l,- 

2  2 

1  1 

—  +  —  i 
2  2 

-1 


(2.95) 


By  inspection,  it  is  evident  that  A2 


=  -l ±i,  satisfying  the  relationship 


defined  in  (2.55).  The  FRFs  for  the  forced  system  follow  (2.58)  for  mode  2.  Note 
that  there  is  a  division  by  zero  in  this  expression  for  the  first  mode  because  of  the 
double  root  at  zero  (multiplicity  is  two),  which  (2.58)  cannot  accommodate. 


It  is  evident  that  the  computations  for  the  general  case  of  viscous  damping 
are  rather  complicated.  The  computations  become  much  simpler  for  a  system 
exhibiting  proportional  damping.  Because  the  example  system  studied  here  was 
specifically  chosen  to  be  proportionally  damped,  the  eigenvectors  are  the  same 
as  for  an  undamped  system,  as  shown  previously.  The  damping  matrix  C  can  be 
transformed  to  the  modal  damping  matrix  of  (2.64)  if  the  following  condition  is 
true: 


(m^1k)(m^1c)  =  (m^1c)(m’1k)  . 


(2.96) 


Evaluating  (2.96)  gives 
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(m_1k)(m_1c) 


(M“1C)(M‘1k) 


1 

0 

-1 

1 

-1 

\ 

1 

0 

-1 

1 

-1 

\ 

2 

-2 

0 

1 

-1 

1 

) 

0 

VL 
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-1 

1 

/ 

-2 

2 

1 

0 

-1 

1 

-1 

\ 
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0 

-1 

1 

-1 

\ 

2 

-2 

0 

1 

-1 

1 

J 

0 

1 

-1 

1 

/ 

-2 

2 

.(2.97) 


Therefore,  the  modal  damping  matrix  is 


C  = 


0  0 
0  4 


So,  for  the  unforced  system, 


2 

0 

q  + 

0 

0 

q+ 

0 

0 

0 

2 

0 

4 

0 

4 

q  =  0 


(2.98) 


(2.99) 


giving 


<7i  =  0 

and 


q2  +  2q2  +  2q2  0  . 


(2.100) 


(2.101) 


The  eigenvalues  are  computed  from  the  characteristic  polynomial.  For  the  first 
mode,  by  inspection  (ol  \  =  0  (shown  earlier).  For  the  second  mode, 


co,  =  co1 


Vl--CF  =  V2  1- 


V2 

v  2  y 


(2.102) 


These  values  match  those  from  the  general  case  computed  earlier.  However, 
(2.102)  shows  that  co2  «.co2  due  to  a  high  modal  damping  value  relative  to  modal 
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stiffness.  As  discussed,  the  mode  shape  matrix  V  is  the  same  as  for  the 
undamped  system  in  (2.77). 

For  comparison,  the  FRF  for  the  first  mass  due  to  the  second  force  is 
constructed  below. 


V  V  V  V 

Hp(ft>)  =  ~  ,  1!  21  2_  +-  .  !:  22  2  _ 

Kj  +  icocl  -  co  ml  k2  +  icoc2  -  co  m2 

__L _ 1 

-2  co2  A  +  Aico-lco2 


(2.103) 


The  impulse  response  to  the  first  mode,  by  inspection,  is  linear  towards  infinity, 
so  it  provides  no  new  insight.  The  presence  of  damping  greatly  changes  the 
behavior  of  the  second  mode  (compare  Figure  7  with  Figure  6(c)).  The 
oscillations  damp  out  quickly  because  the  damping  is  large.  Notice  in  Figure  7 
how  the  damping  roughly  halves  the  peak  amplitude  of  oscillation  of  the 
undamped  system. 


Figure  7.  Proportionally  damped  impulse  response  of  2H12  for  c  =  1 , 

c2  =  4. 
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As  an  illustration  of  the  effect  of  the  relationship  between  the  damping  and 
stiffness  constants,  Figure  8  shows  the  impulse  response  of  the  second  mode 
when  the  value  of  damping  is  decreased  by  an  order  of  magnitude  to  c  =  0.1, 
while  keeping  k  =  1.  Less  damping  increases  the  settling  time  of  the  system  while 
also  increasing  its  peak  amplitude. 


Figure  8.  Proportionally  damped  impulse  response  of  2H12  for  c  =  0.1, 

c2  =0.4. 

G.  FORCED  VIBRATION  WITH  VISCOUS  DAMPING  IN  THE  LAPLACE 
DOMAIN 

Modal  analysis  can  also  be  completed  in  the  Laplace  domain,  where 
manipulation  of  the  equations  is  easier.  Many  conventional  control  techniques 
have  deep  roots  in  the  Laplace  domain.  To  work  in  the  Laplace  domain,  simply 
substitute  co  =  into  the  FRFs  discussed  previously. 

Assuming  zero  initial  conditions,  Equation  (2.1)  describing  the  motion  of 
the  system  in  Figure  3  becomes  Z(s)X(s)  =  F(s),  where  Z(s)  =  Mr  +Cs  +  K . 

Matrix  ZeM“"  is  called  the  dynamic  stiffness  matrix  [9],  The  transfer  function 
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HgR"m  between  the  displacement  XeM"xl  and  applied  force  FeMnxl  is  the 
inverse  of  the  dynamic  stiffness  matrix  and  is  given  by 


H(s)  =  [Ms2  +  Cs  +  K]~‘ 


N (5) 

d(s ) 


adj[Mr +Cs  +  K] 
det[Ms2  +  Cs  +  K]  ' 


(2.104) 


In  (2.104),  d(s )  is  the  characteristic  polynomial,  whose  roots  are  the 

complex  conjugate  pole  pairs,  of  the  system,  and  NeK"x"  is  a  matrix  of 
numerator  polynomials  that  relates  to  residues  (discussed  below).  The  roots  of 
the  characteristic  polynomial  are  the  same  as  the  eigenvalues  of  the  chain  of 
single-DOF  systems,  discussed  previously.  Recall  from  (2.68)  that  Xr  contains 
both  a  real  part  and  imaginary  part,  such  that  Xr  =  o)£r  ±  ico'r .  The  roots  of  the 
numerator  polynomials  are  the  zeros  of  the  system  [1 3], 

The  transfer  function  in  (2.104)  can  be  deconstructed  into  components 
separately  describing  motion  of  the  system  as  a  whole  (termed  “rigid-body”)  and 
the  oscillatory  motion  (“flexible-body”)  using  partial  fraction  expansion  but  only  if 
the  system,  indeed,  exhibits  rigid-body  motion.  A  more  general  case  is  discussed 
later.  These  components  represent  the  modes  from  the  previous  section.  For  the 
2-DOF  system  shown  in  Figure  3,  we  have  (from  (2.104)) 


H(s)  = 


m2s~  +cs  +  k 
cs  +  k 


cs  +  k 


m{s  +cs  +  k 


s2  {jnlm2s>  +  (m,  +  m2 )  cs  +  (m,  +  m2)kj 


(2.105) 


Using  partial  fraction  expansion  on  each  element  of  H(s)  allows  the 
individual  transfer  functions  between  the  bodies  to  be  determined.  Let  H.,(.s-) 
equal  the  transfer  function  for  the  displacement  Xjk(s)  and  FA  (s),  such  that 
Xjk(s)  =  H^XsjFjXs) .  In  other  words,  Xn(,!>)  is  the  displacement  of  Xj  due  to  /j 
and  Xl2(s)  is  the  displacement  of  x,  due  to  f2 .  From  (2.105), 
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H  „(*)  =  -  + 


s2  mlm2s2  +  (m,  +  m2)cs  +  (mi  +  m2)k 

The  numerator  of  (2.106)  is 

a[my  m2s2  +  (m,  +  m2)cs  +  (m[  +m2)ksj  +  bs 2  =  m2s2  +cs  +  k 

Solving  for  a  and  b  we  obtain 


(2.106) 


(2.107) 


a  = - 

777,  +  777, 

,  777,777, 

0  =  777, - 1 - — 

777,  +  777, 


(2.108) 


Substituting  (2.108)  back  into  (2.106)  gives 


777, 


H„(s)  = 


■  + 


777,  +  777, 


(  777,  +  777,  )  52  777, 777,  5 ?  +  ( 777 ,  +  777,  )  C.V  +  ( 777 ,  +  777,  )  k 


(2.109) 


The  first  term  in  (2.109)  represents  the  rigid-body  motion  of  mass  1  (mode  1) 
while  the  second  term  represents  the  flexible-body  motion  (mode  2)  due  to  its 
attachment  to  mass  2.  Similarly,  the  remaining  transfer  functions  are  given  by 


H„(s) 


1 

(  777,  +  777,  )  .S' 2 


777,777, 

777,  +  777, 

777|777,52  +  ( 777,  +  777,  ) C.S  +  (t77,  +  777,  ) k 


H,,(s) 


H,,(s) 


777,777, 

1  777,  +  777, 

(?77|  +  777,  ) .S'2  777 1 777, S’ 2  +  (777,  +  777,  ) C.V  +  (777,  +  777,  ) /c 


1  777,  +  777, 

(777,  +  777, ) .S'2  777, 777, 5 2  +  (/77,  +  777,  ) CS'  +  (t77,  +  777,)/: 


(2.110) 


The  symmetry  is  expected,  as  the  motion  of  the  two  masses  occurs  only 
along  the  same  axis,  and  they  are  only  connected  to  each  other.  The  first  terms 
in  each  transfer  function  correspond  to  the  rigid-body  motion  of  the  system.  They 
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contribute  to  the  first  mode  where  co{  =0 .  The  second  terms  model  flexible-body 
motion  and  represent  the  second  mode. 

The  residues  of  these  functions  are  proportional  to  the  modal  vectors  and 
can  be  used  to  construct  the  residue  matrix  RreC“  [11].  (“Residue”  is 
analogous  to  “modal  constant”  shown  previously.)  In  general,  assuming  all  poles 
have  multiplicity  of  one  [9], 


H(s)  =  £ 


r=l 


R 


R 


S-J 5  -  A, 


V  J 


(2.111) 


where  A,.  and  X*  are  the  complex  conjugate  pole  pairs  and  n  is  the  number  of 
modes  (or  degrees  of  freedom).  Residues  are  computed  using  partial  fraction 
expansion,  as  above.  The  mode  shape  of  mode  r  is  described  by  a  matrix 

'FeC"M  whose  columns  are  normalized  modal  vectors,  y/r,  proportional  to  the 
corresponding  elements  of  R,. .  As  was  shown  in  the  case  of  a  proportionally 
damped  system,  the  modal  vector  y/r  equals  vr  from  the  undamped  case.  The 
modal  vectors  are  calculated  by  using  the  residues  from  one  row  or  column  of  H 
(for  both  modes).  Interestingly,  since  the  system  is  proportionally  damped,  the 
residues  will  be  purely  imaginary  [1 1], 

The  inverse  Laplace  transform  of  (2.111)  gives  the  impulse  response 
function  h (t)  studied  earlier. 


hM =£(]*/■ 

r= 1 


(2.112) 


H.  NUMERICAL  EXAMPLE  IN  THE  LAPLACE  DOMAIN 


Assume  a  nominal  case  where  mt  =  m2  =  c  =  k  =  \ .  For  the  undamped 


case,  Aj  =  0 ,  A2  =  2  ,  and  V  = 


-1 

1 


(from  (2.10)  and  (2.16)),  and  it  follows 


that  the  natural  frequencies  are  co{=  0  and  =  Eigenvector  v2  indicates 
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that  the  amplitudes  of  motion  of  both  masses  in  mode  2  will  be  equal  but  in 
opposite  directions. 


The  equations  of  motion  in  the  Laplace  domain  become 

f  1  1 


X»  = 
X12(s)  = 
X21(s)  = 
x22(5)  = 


V5 

n 

2 


VS 


2  +  s2+2s  +  2 


1 


s  +  2s  +  2  J 


Vf 

J  2  F| 

'IF 

2  2 


1 


s  +  2s  +  2 


2Fl 


,  2  2 
V5 


1 


s  +  2  s  +  2  / 


—  Ft 


(2.113) 


It  is  readily  apparent  that  the  poles  of  the  rigid-body  component  (r  =  l)  are 
at  s  =  0.  The  rigid-body  mode  from  the  double-integrator  terms  in  each  equation 
implies  that  the  system  moves  as  a  whole  in  the  same  direction  due  to  the 
applied  forces.  The  system’s  second  mode  comes  from  the  oscillatory  motion 
between  the  two  masses,  modeled  in  the  second  terms.  Take  the  transfer 
function  component  of  the  second  term  of  X12(s).  The  poles  and  corresponding 

residues  are 


X2  =  -\  +  i 
2R12  =  ±  — 


(2.114) 


Recall  that  the  poles  match  the  eigenvalues  calculated  earlier.  The  residues  can 
be  used  to  reconstruct  the  transfer  function  from  (2.111 ). 


Hi2(5) 


1  . 
—i 

4 


s  -(-!  +  /) 


1  . 
—i 

4 


s-(-l-i) 


2 

s~  +  2s  +  2 


(2.115) 


The  residue  matrix  may  be  completed  by  following  the  same  procedure. 
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(2.116) 


R, 


_ 1  .  ,1. 

+— i  ±—i 

4  4 

,  1  •  _1. 

±  —  i  +  —  i 

4  4 


The  corresponding  impulse  response  function  is 


2h12(f)  =  —  ie'  ie(  1  '^  =  — —  e  'sinf  . 

2  12  4  4  2 


(2.117) 


The  time  history  of  (2.1 17)  is  precisely  the  same  as  the  plot  shown  in  Figure  7. 


I.  MASS-SPRING-DAMPER  SYSTEM  FIXED  TO  A  STRUCTURE 

Suppose  now  that  the  system  from  Figure  3  is  connected  on  one  side  to  a 
fixed  structure,  as  in  Figure  9.  This  connection  represents  a  controller  that  may 
be  implemented  to  stabilize  the  rigid  motion  of  the  flexible  system. 


;  m 


f2(t) 


=  0 


x2  =  0 


Figure  9.  Two-DOF  mass-spring-damper  system. 


The  first  mode  of  the  system  no  longer  exhibits  rigid-body  motion  since  the 
masses  are  fixed  in  place  by  the  first  spring-damper  pair.  The  new  equations  of 
motion  are 


/j  =  mx xx  +  (c,  +  c2 )i,  +  (&,  +k2)xx- cx2  - k2x2 
f2  =  m2x2  -  c2xx  -  k2xx  +  c2x2  +  k2x2 

which  leads  to 


(2.118) 
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.  (2.119) 


Mx(0  +  C  x(t)  +  K x(t)  =  f  (t) 


r  m,  0 

C|  +  c2  —c0 

ky  +  k2  —k2 

o 

to 

x(l)  + 

x(l)  + 

x(r)  =  f(0 

1 

CN 

o 

(N 

1 

_ 1 

~k2  k2 

Assuming  the  same  nominal  values,  ml  =  m2  =cl=c2  =  kl=k2=  1 ,  the 
natural  frequencies  and  mode  shapes  for  the  undamped,  unforced  system  are 
given  as 


[Ar]  =  [©/]  = 


0.382  0 

0  2.618 


V  = 


0.618 

1 


-1.618 

1 


(2.120) 


From  (2.120),  the  modal  matrices  are  obtained  as 


0.618  1 

1 

0 

0.618  - 

-1.618 

1.382 

0 

-1.618  1 

_  0 

1 

1 

1 

0 

3.618 

0.618  1 

2 

-1 

0.618 

-1.618 

0.528 

0 

-1.618  1 

-1 

1 

1 

1 

0 

9.472 

0.618  1 

2 

-1 

0.618 

-1.618 

0.528 

0 

-1.618  1 

-1 

1 

1 

1 

0 

9.472 

(2.121) 


The  equation  of  motion  in  modal  coordinates  is,  therefore, 


1.382 

0 

q+ 

0.528 

0 

q+ 

0.528 

0 

q  = 

0.618  1 

0 

3.618 

0 

9.472 

0 

9.472 

1 

t—H 

00 

t—H 

VO 

H 

1 

_ 1 

The  equations  for  each  mode,  after  dividing  through  by  mr ,  are  thus 

<7,  +  0.382^  +  0.382^  =0.447/,  +  0.724/2  =  (2.122) 

1 .3o  Z 

q2  +  2.6184,  +  2.618^  =-0.447/,  +0.276/,  =  -^~  .  (2.123) 

3.618 


These  equations  describe  the  motion  of  two  single-DOF  systems,  where  the 
motions  of  the  modal  masses  are  uncoupled  and  are  of  the  form 
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qr  +  2C,.COrqr  +  COr2q , 


(2.124) 


For  the  forced  system, 


col  =  V0.382  =0.618 
0.382 


c 


(2)(0.618) 


0.309 


co'  =  0.588 


and 


(02 

c2 


=  a/2.618  =  1.618 
2.618 


(2)(1.618) 


=  0.809 


ffl'=0.951 


(2.125) 


(2.126) 


Equations  (2.125)  and  (2.126)  can  be  interpreted  as  follows:  mode  1,  with 
the  smaller  natural  frequency,  represents  the  main  motion  of  the  system, 
whereas  mode  2,  with  the  larger  natural  frequency,  represents  the  flexible  effect. 
This  flexible  effect  perturbs  the  main  motion  of  the  system.  In  this  context,  it  is 
desirable  to  minimize  the  magnitude  of  q2  and  its  derivatives  since  the  flexibility 
is  parasitic  to  the  performance  of  the  control  system.  This  result  illustrates  the 
primary  reason  why  modal  analysis  is  a  valuable  tool  in  control  system  design. 

The  controlled  2-DOF  system  shown  in  Figure  9  cannot  be  deconstructed 
into  rigid-body  and  flexible-body  components,  as  shown  in  the  previous  section, 
since  motion  of  both  masses  is  oscillatory.  Assuming  the  same  nominal  values 
m,  =  m2  =  c,  =  c2  =  kx  =  k2  =  1 ,  in  the  Laplace  domain  the  FRF  is 


s2  +5  +  1 

5  +  1 

N(s)  adj[M,r  +CX  +  K.] 

5  +  1 

5’  +  25  +  2 

~  d(s)  ~  det[Mr  h-C^  +  k]  ' 

54  +  353+452+25  +  l 

(2.127) 


where  the  poles  are  \  =  -0.191  ±0.588/  and  A2  =-1.311  ±0.95 li ,  which  match 
the  values  computed  previously.  The  residue  matrices  for  each  mode  are 
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(2.128) 


+0.235/  +0.380/ 
+0.380/  +0.615/ 

+0.380/  ±0.235/ 
±0.235/  +0.145/ 


Therefore,  the  mode  shape  matrix,  using  the  first  row  of  each  residue  matrix,  is 


R11 

2R11 

R2i 

2^21 

R2i 

D 

2^21 

R2i 

D 

2^21 

-0.235  -0.380 

-0.380  0.235 

1  1 


0.618  -1.618 

1  1 


(2.129) 


While  determining  mode  shapes  and  residues  in  the  Laplace  domain  is 
straightforward,  the  transfer  functions  in  H  cannot  be  used  to  get  the  modal 
coordinates  of  the  system.  They  do,  however,  provide  meaningful 
representations  of  the  behavior  of  the  system  in  physical  space. 


J.  SUMMARY 

The  motion  of  flexible  multi-body  systems  is  complicated  due  to  the 
interaction  between  the  various  degrees  of  freedom.  Modal  analysis  enables  the 
description  of  the  motion  to  be  decomposed  into  a  series  of  single-DOF  systems, 
each  with  different  modal  properties  (i.e.,  frequency  and  damping).  Control 
applications  are  typically  concerned  with  manipulating  the  behavior  of  the 
dominant  (slow)  modes  while  minimizing  the  excitation  of  the  flexible  (fast) 
modes.  Thus,  modal  coordinates  are  a  natural  space  for  synthesis  of  rapid  slew 
maneuvers,  and  this  chapter  has  provided  the  necessary  background  for 
accomplishing  this  task. 
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III.  CONVENTIONAL  CONTROL  OF  A  3-DOF  SYSTEM 


This  chapter  increases  the  complexity  of  the  mass-spring-damper  system 
by  adding  a  third  DOF.  The  system  is  now  analogous  to  a  controlled  spacecraft 
body  with  a  flexible  antenna  mounted  on  a  double-axis  gimbal.  The  controller  is 
used  to  maintain  the  spacecraft’s  position  and  move  the  links  of  the  notional 
antenna.  This  chapter  analyzes  a  model  to  understand  the  behavior  of  such  a 
system  from  the  perspective  of  conventional  control  techniques. 

A.  DEFINING  THE  MODEL 

Consider  the  three-mass  system  shown  in  Figure  10.  It  is  representative  of 
a  controlled  spacecraft  with  a  2-DOF  flexible  appendage.  In  this  case,  the 
stiffness  k{  and  damping  c,  can  be  attributed  to  the  feedback  control  of  the  base 
body,  whereas  k2 ,  k3 ,  c2 ,  and  c3  are  the  stiffness  and  damping  constants  of  the 
flexible  appendage  (locked  in  a  specific  orientation). 


x-t  =  o  x2  =  o  x3  =  o 


Figure  1 0.  Mass-spring-damper  system  with  three  degrees  of  freedom. 

The  equations  of  motion  governing  this  system  are 

/  =  m1x[  +  (q  +  c2 )i,  —  c2x2  +  [kx  +  k2)x ,  -  k2x2 

f2  =  m2x2  -  c2x j  +  (c2  +  c3 )  x2  -  c3i3  -  k2x]  +  ( k2  +  k3 )x2-  k3x3  .  (3 . 1 ) 

f3  =  m3x3  -  c3x2  +  c3x3  -  k3x2  +  k3x3 
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Thus, 


ml 

0 

0 

c,+c2 

-c2 

0 

0 

m2 

0 

c  = 

-c2 

C2  +C3 

-c3 

0 

0 

m3 

0 

-c3 

c3 

k j  +  kj 

—k2 

0 

Xt(t) 

"  f(0 

-k2 

k2  +  k3 

-k3 

-  x(/)  = 

x2(t) 

,  and  f(/)  = 

fit) 

.  Unit  values  for 

0 

^3 

x3(t) 

1 

1 _ 

each  parameter  in  the  system  are  assumed  for  convenience  so  that 

ml  =  m2  =  m3  =  c,  =  c2  =  c3  =  k{  =  k2=k3  =  1 . 


B.  ANALYSIS  OF  THE  COUPLED  SYSTEM 

Using  the  method  described  in  Chapter  II,  the  transfer  function  matrix  for 
the  system  in  Figure  10  can  be  written  as 


H  (s) 


s4  +  3  s  +  4  s~  +  2  s  +  1  s  +  2  s~  +  2s  +  1  s~  +  2s  + 1 

s3  +  2s2  +  2s  +  1  s4  +  3s3  +  5s2  +  4s  +  2  s3  +  3s2  +  4s  +  2 

s2  +  2s  + 1  s3  +  3s2  +  4s  +  2  s4+4s3  +  7s2  +  6s  +  3 

s6  +  5s5  +  11s4  +  13s3  +  9s2  +  3s  +  1 


(3.2) 


Recall  that  the  poles  of  this  transfer  function  matrix  are  the  roots  of  the 
denominator  polynomial  (or  “characteristic  equation”),  which  are  equal  to  the 
eigenvalues  of  the  system 

[A,.]  =  diag[  -0.990 ±0.4339/  -0.7775  ±0.9749/  -1.6235  ±0.7818/  ].  (3.3) 


Using  the  procedure  discussed  in  Chapter  II,  the  undamped  mode  shape  vectors, 
rounded  to  four  decimal  places,  are 


1  1  1 
1.8019  0.4450  -1.2470 

2.2470  -0.8019  0.5550 


(3.4) 


The  modal  matrices,  then,  are 
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(3.5) 


M  =  VrMV  =  diag[  9.2959  1.8412  2.8629  | 

C  =  VrCV  =  diag[  1.8412  2.8629  9.2959  ]  . 

K  =  VrKV  =  diag[  1.8412  2.8629  9.2959  j 

From  these  equations,  the  natural  and  damped  frequencies,  cor  and  co'  (in 
radians  per  second),  as  well  as  the  critical  damping  ratios,  ,  may  be  computed 
as 

[«,.]  =  diag[  0.4450  1.2470  1.8019  j 

[£,.]  =  diag[  0.2225  0.9749  0.7818  ]  .  (3.6) 

[®,']  =  diag[  0.4339  0.9749  0.7818  j 

The  relative  values  of  the  modal  masses  in  (3.5)  help  to  illustrate  to  an 
important  property  of  this  system.  Modal  mass  m,  =  9.2959  is  five  times  larger 
than  rh2  and  over  three  times  larger  than  m3.  The  damping  ratios  increase 
through  the  modes.  One  could  hypothesize,  then,  that  the  motion  of  mass  1 
dominates  the  total  motion  of  the  system,  just  as  was  observed  in  the  system 
with  rigid-body  motion  (see  Chapter  II).  Analysis  of  the  FRFs  and  IRFs  will 
support  this  hypothesis.  Figure  1 1  shows  the  Bode  frequency  response  and 

Figure  12  shows  the  motion  of  mass  1  in  response  to  impulse  force  1  for  all 

modes. 
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Figure  1 1 .  Bode  magnitude  and  phase  plot  of  Hn  (to) 

for  a  3-DOF  system. 


In  Figure  1 1,  the  vertical  dashed  and  dot-dashed  lines  represent  the  damped  and 
natural  frequencies  for  each  mode,  respectively.  Clearly,  the  total  motion  of  mass 
1  due  to  force  1  is  greatly  impacted  by  mode  1,  especially  in  the  low  frequency 
range.  This  is  because  the  magnitudes  of  mode  2  and  mode  3  are 
-10  dB  and  -20  dB,  respectively,  in  the  low  frequency  range.  In  the  high 
frequency  range,  the  responses  of  modes  2  and  3  dominate.  This  is  why  flexible 
systems  oscillate;  sudden  changes  in  the  applied  force  excite  these  high 
frequency  effects.  The  impulse  responses  (Figure  12)  confirm  that  mode  1 
dominates  the  response  and  that  modes  2  and  3  decay  quickly  due  to  the  large 
damping.  If  this  damping  is  reduced,  the  higher  modes  will  be  superposed  onto 
mode  1  to  a  greater  degree.  This  is  the  characteristic  that  degrades  the 
performance  of  a  control  system  and  should  be  avoided  via  a  combination  of 
control  and  reference  trajectory  design.  Figure  13  illustrates  this  point;  it  shows 
the  IRF  hn (r)  for  reduced  damping  (c,  =  c2  =  c3  =  O.l) . 
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Figure  12.  Impulse  response  function  hn(r)  for  a  3-DOF  system  with 

“heavy”  damping  (cj  =  c2  =  c3  =  l) . 


Figure  13.  Impulse  response  function  hu(r)  for  a  3-DOF  system  with 

“light”  damping  (c,  =  c2  =  c3  =  O.l) . 
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Figure  14  shows  the  impulse  responses  for  the  transfer  function  matrix  in 
(3.2).  The  repeated  impulse  responses  (as  for  H12  and  H2I,  for  example)  result 
from  the  symmetry  of  the  matrix  H(s)  .  It  is  apparent  from  the  figure  that  m2  and 
m3  experience  the  greatest  displacement  about  their  equilibrium  positions  due  to 
the  forces  acting  on  them.  This  follows  due  to  the  nature  of  the  system;  each 
mass  is  less  constrained  than  its  previous  counterpart  as  one  moves  away  from 
the  fixed  structure  (base).  In  the  control  of  a  multi-body  system,  this  is  important 
because  if  forces  2  and  3  are  not  carefully  controlled,  they  can  cause  large 
displacement  of  the  base  body  (i.e. ,  mass  1).  See  responses  for  Hl2  and  Hl3  in 
Figure  14. 


Figure  14.  Impulse  responses,  h/7.(r),  for  system  in  Figure  10. 
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The  residue  matrices  for  each  mode,  using  partial  fraction  expansion  from 
Chapter  II,  are 


+0.1240/ 

+0.2234/ 

+0.2786/ 

+0.2234/ 

+0.4025/ 

+0.5019/ 

+0.2786/ 

+0.5019/ 

+0.6259/ 

+0.2786/ 

+0.1240/ 

±0.2234/ 

+0.1240/ 

+0.0552/ 

±0.0994/ 

±0.2234/ 

±0.0994/ 

+0.1791/ 

+0.2234/ 

±0.2786/ 

+0.1240/ 

±0.2786/ 

+0.3473/ 

±0.1546/ 

+0.1240/ 

±0.1546/ 

+0.0688/ 

(3.7) 


These  residues  may  be  used  to  construct  mode  shape  vectors 


iRn 

2  R1 1 

3R1 1 

iRn 

2^11 

3^11 

1R21 

D 

2^21 

D 

3^21 

1R11 

2  R1 1 

3^11 

1  R31 

2  R31 

p 

3^31 

1R11 

2  R1 1 

3R1 1 

1  1  1 
-0.2234/  -0.1240/  0.2786/ 

-0.1240/  -0.2786/  -0.2234/ 

-0.2786/  0.2234/  -0.1240/ 

-0.1240/  -0.2786/  -0.2234/ 


1  1  1 
1.8019  0.4450  -1.2470 

2.2470  -0.8019  0.5550 


(3.8) 


Figure  15  plots  these  mode  shapes.  Recall  that  a  plot  like  Figure  15  is  useful  in 
understanding  the  relative  motions  of  each  degree  of  freedom  of  the  system  for 
each  modal  frequency.  Each  line  in  Figure  15  represents  a  single  mode  mapped 
to  each  DOF.  For  instance,  motion  of  mode  1  (blue  line)  causes  all  three  masses 
to  move  in  the  same  direction  (because  each  amplitude  has  the  same  sign)  but 
with  varying  magnitudes.  One  should  expect  m3  to  oscillate  at  more  than  double 
the  amplitude  of  m1  (2.25:1  ratio)  and  only  slightly  more  than  m2  (1.80:1  ratio). 
Mode  2  (red  line)  would  see  ra3  moving  in  a  direction  opposite  the  other  two 
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masses  because  its  amplitude  has  the  opposite  sign  of  the  other  two,  and  so  on. 
(It  can  also  be  shown  that  these  mode  shapes  can  be  formed  from  the  peak 
amplitudes  of  the  imaginary  parts  of  their  FRFs.) 


Figure  1 5.  Mode  shapes  of  3-DOF  system  (c,  =  c2  =  c3  =  l) . 

C.  MODELING  A  SPACECRAFT  SLEW 


With  minor  modifications,  Figure  10  can  be  used  as  a  representative 
model  for  a  spacecraft  antenna  slew.  Suppose  a  steerable  antenna  is  attached  to 
the  body  of  a  spacecraft  (see  Figure  16)  by  a  two-axis  gimbal.  For  simplicity, 
assume  the  center  of  rotation  of  the  gimbal  is  aligned  with  its  center  of  mass,  and 
the  antenna  is  mounted  in  such  a  way  that  there  is  no  offset  between  it  and  the 
spacecraft  body. 
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Figure  16.  Notional  schematic  for  a  two-axis  antenna  mounted  on  a 

spacecraft  body,  after  [5], 

The  system  in  Figure  10  in  its  current  form  can  be  used  as  a  notional 
model  for  the  flexible  component  of  motion  (vibrations)  resulting  from  the  forces 
applied  to  the  azimuth  and  elevation  bodies,  where  the  spacecraft  body  is 
analogous  to  mx,  azimuth  to  m2 ,  and  elevation  to  m3.  The  spring  and  damper 
connecting  mass  1  to  the  rigid  structure  represents  the  spacecraft  controller  that 
seeks  to  counteract  the  motion  induced  by  the  antenna  and  maintain  the 
spacecraft  in  a  certain  position.  The  forces  applied  to  the  antenna  to  complete 
the  slew  are  the  forces  that  are  input  into  the  equations  of  motion  in  (3.1).  The 
forces  are  also  the  inputs  to  simulations  of  the  transfer  function  H(s)  in  (3.2). 

In  order  to  model  the  entire  maneuver — not  only  the  flexible  motion  in  the 
notional  antenna’s  locked  position — Figure  10  must  be  modified  to  include  the 
gross  motion  of  the  notional  gimbals.  A  revised  model  is  shown  in  Figure  17.  The 
model  shows  that  masses  2  and  3  have  moved  to  new  equilibrium  positions 
defined  by  the  desired  final  position  at  time  tf ,  while  the  flexible-body  motion  has 

been  preserved  and  referenced  to  a  non-stationary  equilibrium  point. 
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*l(ff)  ^  X2(ff)  -  X2f  *3(ff)  *3f 

(b) 

Figure  17.  Three-DOF  system  modeling  rigid  translation  with  flexible 


effects:  (a)  initial  position;  (b)  final  position. 

Let  X2  =  x2(t)-x2{t0)  and  X3  =  x3(t)- x3(t{))  be  the  translation  of  each  mass  as  the 
prismatic  joints  are  moved.  The  motion  of  the  masses  now  also  includes  what 
can  be  equated  to  rigid-body  motion.  The  total  motion  of  the  masses  is  the  sum 
of  this  rigid-body  motion  and  the  flexible-body  component  for  each  mass.  Thus, 
the  equations  for  the  total  motion  become 

/,  =  mxX j  +  m,  x,  +  (c,  +  c2)xl- c2x2  +  (k]  +  k2 ) x]  - k2x2 

f2  =  m2X2  +  m2x2  -  c2i,  +  (c2  +  t’3)i2  -  c3x3  -  k2x{  +  (k2  +  k3)x2  -  k3x3  (3.9) 
f3  —  m3X3  +  m3x3  —  c3x2  +  c3x3  —  k3x2  +  k3x3 

where  x,  =Xx  =X,  =0  since  the  translation  of  the  base  is  zero.  The  total  motion 
of  the  system  is  the  linear  superposition  of  the  rigid-body  motion  and  the 
vibrational  effects. 
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D. 


SLEWING  WITH  COMPUTED  TORQUE  CONTROL 


Assume  the  following  boundary  conditions  and  constraints  for  a  notional 
antenna  slew  are: 

(x1(r0),x2(r0),X3(r0),x1(r0),x2(r0),x3(r0))  =  (o,-lO,5,0,0,0.l) 
(x1(r/),x2(r/),x3(r/),x1(r/),x2(r/),x3(r/))  =  (0,l5,-l5,0,0.l,0) 

M  =  0.75  ,  (3.10) 

|/J  =  l-5 

(^0,^)  =  (0,38.1) 

where  X2  and  X2  are  the  translation  and  velocity,  respectively,  in  the  notional 

azimuth  plane,  and  X3  and  X2  are  the  translation  and  velocity,  respectively,  in 
the  notional  elevation  plane.  Both  planes  are  fixed  relative  to  the  spacecraft  body 
frame.  The  parameters  in  Figure  17  are,  again,  assigned  unit  values. 

Common  practice  for  implementing  a  slew  maneuver  to  satisfy  (3.10)  is  to 
use  a  so-called  bang-bang  or  bang-off-bang  trajectory  where  the  maximum 
allowable  force  is  applied  over  a  short  period  of  time  to  accelerate  each  movable 
joint,  then  the  same  magnitude  of  force  is  applied  in  the  opposite  direction  to 
decelerate  the  links  to  their  desired  positions.  This  approach  leads  to  the  idea  of 
computed  torque  control,  in  which  required  forces  are  generated  from  knowledge 
of  the  current  state  and  its  comparison  to  a  desired  profile.  The  control  system 
uses  feedback  about  the  current  state  to  determine  the  amount  of  control  to 
apply  [4],  By  this  method,  the  spacecraft  can  maneuver  along  the  prescribed 
path.  A  control  law  for  such  a  system  takes  the  form  given  in  (3.11),  which  is  a 
simple  PD  controller. 

f  =  KP{xd~~x)  +  Kd(xdmx)  ■  (3-11) 

In  (3.11),  Kp  and  Kd  are  gains  that  can  be  programmed  to  achieve  desired 

characteristics  of  the  system.  This  method  most  likely  produces  a  suboptimal 
solution  but  is  easy  to  implement.  Figure  18  shows  the  results  of  the  maneuver 
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implementation  using  computed  torque  control  for  a  rigid-body  model;  it  also 
shows  the  behavior  of  the  system  after  completion  of  the  slew  (indicated  by  the 
circles).  Gains  for  the  system  are  based  on  a  one-second  settling  time  with 
critical  damping  «  =  1)  to  minimize  overshoot.  The  control  law  in  (3.1 1)  is  applied 
to  both  axes  (azimuth  and  elevation)  independently  to  achieve  the  desired  final 
conditions  in  (3.10)  within  the  desired  time  of  38.1  seconds. 
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Figure  18.  Rigid-body  slew  maneuver  using  computed  torque  control:  (a) 

position;  (b)  rate;  (c)  acceleration. 

To  model  the  flexible-body  motion  of  the  system,  the  responses  of  the 
transfer  function  matrix  in  (3.2)  can  be  simulated  by  applying  the  forces  produced 
from  the  control  law  in  (3. II).2  Recall  that  X(s)  =  H(s)F(s),  or 

X.(s)  =  H.,(5)Ffc(s)  for  j  =  1,2,3  and  k  =  1,2,3  .  Therefore, 


2  MATLAB's  “Isim”  command  is  very  useful  for  this. 
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(3.12) 


JCl(5)  =  H11(5)FI(5)  +  H12(5)F2(s)  +  H1,(5)F3(s) 

X2(s)  =  H2i(s)Fi(s)  +  H22(5)F2(5)  +  H23(s)F3(s) 

x3(s)  =  'R31(s)F1(S)  +  H32(s)F2(s)  +  H33(s)F3(s) 

We  are  only  concerned  at  this  point  with  x2  and  x3  (the  control  law  will  be  set  to 
maintain  x1  =0).  Figure  19  shows  the  results. 


Figure  19.  Flexible  motion  due  to  slew  in  Figure  18. 


Note  how  the  peak  amplitudes  of  the  vibrational  components  coincide  with  the 
changes  in  acceleration  in  Figure  18,  though  the  vibrations  occur  on  a  much 
smaller  scale  than  the  gross  motion  of  the  system.  In  addition,  there  is  some 
flexible  spillover.  In  other  words,  even  though  the  gross  motion  is  fixed  beyond 
tf  =  38.1 ,  the  system  continues  to  oscillate  in  accordance  with  (3.12). 

Suppose  now  that  it  is  desired  to  constrain  vibrations  to  within  0.1  position 
units  for  both  links.  A  difficulty  in  using  computed  torque  control  arises  because 
one  must  set  the  appropriate  limits  on  force  and  velocity  to  ensure  this  constraint 
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can  be  met.  Finding  these  limits  requires  “tweaking”  the  specifications  of  the 
system  until  the  desired  effect  is  achieved.  However,  the  total  slew  time  can 
increase,  which  may  not  be  desirable.  For  example,  an  increase  in  slew  time  of 
three  seconds  to  41.1  seconds  reduces  vibrations  to  below  the  constraint  (see 
Figure  20).  Figure  21  compares  the  acceleration  commands  for  each  case;  the 
controller  applies  less  torque  during  the  slower  maneuver,  and  this  is  the 
mechanism  used  to  reduce  the  magnitudes  of  the  oscillations.  Hence,  using  a 
conventional  computed  torque  strategy,  the  only  way  to  control  flexible  effects  is 
by  manipulating  the  acceleration  of  the  joints.  This  is  why  a  flexible  antenna  must 
be  slewed  slowly  to  avoid  such  disturbances. 


Figure  20.  Reducing  the  flexible  effects  by  increasing  the  slew  time. 
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Figure  21 .  Comparison  of  acceleration  for  different  slew  maneuver  times. 

E.  SUMMARY 

In  this  chapter,  the  mass-spring-damper  model  was  expanded  to  include  a 
third  degree  of  freedom.  This  system  is  analogous  to  a  spacecraft  body  with  an 
antenna  mounted  on  a  double-axis  gimbal.  Conventional  control  techniques  were 
then  used  to  perform  a  notional  antenna  slew  to  illustrate  the  resulting  flexible 
effects.  Finally,  it  was  shown  that  when  using  conventional  control  flexible  effects 
can  only  be  reduced  by  increasing  slew  time.  There  is  no  analytical  process  for 
determining  the  time  increase  that  would  properly  constrain  vibrations,  so  the 
process  can  be  tedious.  This  simple  exercise  illustrates  the  state  of  practice; 
meeting  constraints  on  vibrations  is  handled  by  operating  the  system  more 
conservatively.  However,  this  limits  system  throughput  since  slew  time  is  wasted 
time.  In  the  next  chapter,  optimal  control  is  used  to  develop  slew  maneuvers  that 
alleviate  these  issues. 
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IV.  OPTIMAL  CONTROL  OF  A  3-DOF  SYSTEM 


It  is  desired  to  minimize  slew  time  in  order  to  maximize  throughput.  In 
addition,  vibrational  effects  due  to  force  inputs  should  be  limited  to  their  existing 
specifications.  That  is,  moving  faster  should  not  perturb  the  system  any  more 
than  the  preset  solution.  Rather  than  use  conventional  control  to  minimize  flexible 
effects,  and  potentially  increase  slew  time,  it  is  possible  to  determine  an  alternate 
control  to  minimize  both  flexibility  and  slew  time.  This  can  be  done  by 
manipulating  the  system  through  application  of  force  in  such  a  way  that  the 
flexible  properties  of  the  system  may  be  used  to  an  advantage.  Finding  such  a 
solution  can  be  done  by  formulating  and  solving  an  optimal  control  problem. 

Following  [6],  the  problem  must  be  redefined  to  satisfy  the  necessary 
conditions  to  achieve  optimal  control  of  the  system.  First,  the  dynamics  of  the 
system  must  be  determined  in  state-space  form.  The  correct  definition  of  the 
state  and  control  variables  is  crucial  to  a  successful  optimal  control  solution. 
Once  the  states  and  controls  are  properly  defined,  we  can  proceed  to  redefine 
our  problem  statement. 

To  construct  an  optimal  control  problem,  a  cost  functional  /(x(*),u(*))  to 

be  minimized  should  be  identified.  The  cost,  a  function  of  the  state  vector  x  and 
control  vector  u  in  the  system,  consists  of  an  endpoint  cost  (expressed  in  terms 
of  the  final  state  values)  and/or  a  running  cost  (calculated  for  each  time  step). 
Thus,  the  cost  functional  can  be  expressed  as 

'/(x(-)^u(-))  =  £’(x(r/))  + ,p(x(0,u(0^  ■  (4.1) 

For  minimum-time  solutions,  the  cost  functional  is  simply  J  =  tf ,  but  variations  of 
this  simple  objective  function  are  described  below. 
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The  next  step  is  to  define  any  boundary  conditions  for  the  states,  controls, 
and  time.  Any  constraints  on  states  or  controls,  (e.g.,  |/max|  =  1.5),  can  also  be 
applied.  The  redefined  problem  statement  could  be  written  as 


Minimize 
Subject  to 


x  =  f(x(f), u(f)) 


*i(°)  *2(°) 


t0=0 

e(xM)  =  ° 

xL<xi<xu,  i  =  1,2,... n 
uL  <  u,.  <  \xu,  i  =  1,2,... m 


(4.2) 


To  solve  (4.2),  the  Hamiltonian  must  be  derived.  The  Hamiltonian  is  a  function  of 
the  states,  costates,  and  controls,  where  the  vector  of  costates  is  related  to  the 
dynamics  of  the  system  [6],  The  Hamiltonian  is  defined  as 

//(x,A,u)  =  F(x,u)  +  Arf(x,u)  (4.3) 


where  X  are  the  costates. 


The  necessary  conditions  to  achieve  an  optimal  control  solution  can  now 
be  determined  [6],  The  first  necessary  condition  is  the  Hamiltonian  minimization 
condition,  written  as 


c)H 

3u 


=  0  . 


(4.4) 


The  second  necessary  condition  is  the  satisfaction  of  the  adjoint  equation(s) 


i(t) 


dH 

dx 


(4.5) 


The  transversality  condition  relates  the  final  values  of  the  costates  to  the 
Endpoint  Lagrangian  and  comprises  the  third  necessary  condition,  written  as 


(4.6) 
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where  the  Endpoint  Lagrangian  £’(v,x(r/))  =  £’(x(r/))  +  vre(x(r/))  [6], 

The  MATLAB  software  package  DIDO  automatically  determines  an 
optimal  control  solution  that  satisfies  the  above  conditions  [6],  DIDO  is  used  to 
solve  the  optimal  control  problems  posed  through  the  rest  of  this  thesis. 

A.  STATE-SPACE  FORM  OF  COUPLED  SYSTEM 

To  construct  the  optimal  control  problem,  the  system  in  (3.9)  must  be 
transformed  to  state-space  form.  Assuming  unit  values  for  mass,  stiffness,  and 
damping,  let  y  e  M3*1  be  the  total  displacement  of  each  mass  such  that 


Ti 

Xj  +  X j 

y2 

= 

X2  H-  X2 

y3 

X3  +  X3 

Let  z  e  M6j1  be  a  state  vector  such  that 


z  = 


y 

y 


Xj  +  .Xj 

%2  “1“  X-2 

X3  +  X3 
Xj  +  Xj 

%2  X2 
X3  +  X3 


(4.7) 


(4.8) 


Let  u  e  M3j1  be  the  control  vector  u  = 


fi  fi  /3 


I T 


(i.e.,  the  vector  of  forces 


acting  on  each  mass).  The  states  represent  the  total  motion  described  by  the 
sum  of  the  rigid-body  component  and  the  flexible-body  component,  as  described 
by  (3.9).  The  generic  state-space  representation  for  a  rigid-body  slew  is 


X 

X 


0  1 
0  0 


(4.9) 
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where  X  is  the  displacement,  X  =  v  is  the  velocity,  X  =  a  is  the  acceleration,  m 
is  the  mass,  and  f  is  the  control  force.  This  state-space  form,  adapted  for  the 
3-DOF  system,  is 


*1 

Z2 

*3 

Z2 

*3 


0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 


0  0 


(4.10) 


where  u  =  ^  0  f2  /3  J  since  the  restoring  force  is  due  to  the  control  law  on 

mass  1  comprising  kx  and  ci  (which  will  be  accounted  for  as  part  of  the  flexible 
dynamics).  The  flexible-body  component  in  state-space  form  becomes 


(4.11) 
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Developing  the  state-space  representation  using  (4.8)  as  the  state  vector  is 
grueling.  One  cannot  simply  add  the  rigid-  and  flexible-body  state-space 
equations  together  to  produce  a  form  that  is  conducive  to  modal  analysis.  The 
difficulty  lies  in  that  the  flexible  motion  equations  contain  no  X  terms  or  their  time 
derivatives.  For  example,  the  equation  for  z4  is 

z4  =  xx  +  X,  =  2/j  -  2x1  +  x2  -  2i,  +  x2  .  (4.12) 

Instead,  separate  the  rigid-body  and  flexible-body  components  and  define  a  new 
state  vector.  Let  y  e  R12vl  be  the  state  vector 

1T 

y=  Xt  X2  X3  Xx  X2  X3  a,  a,  a3  a,  a2  i3  (4.13) 

and 

-i  T 

y  =  Xx  X2  X3  X,  X2  X3  ij  x2  x3  xx  x2  x3  .  (4.14) 

The  state-space  representation,  with  unit  values  substituted,  is  therefore 

0001000  0  0  0  0  0  1  [o  0  0 

000010000000  000 

000001000000  000 

000000000000  100 

000000000000  010 

0000000  0  0  0  0  0  001  (4  15) 

0000000  0  0  1  0  0  000 

0000000000  1  0  000 

00000000000  1  000 

000000  -2  1  0  -2  1  0  100 

0000001  -2  1  1  -2  1  010 

0000000  1  -1  0  1  -lJ  |_  0  0  1 

Equation  (4.15)  provides  the  dynamics  required  for  the  optimal  control  problem 
definition.  Recall  that  (4.15)  represents  coupled  motion  in  physical  space.  The 
time-optimal  slew  given  the  conditions  in  (3.10)  and  dynamics  in  (4.15)  can  be 
determined  by  solving  the  following  problem: 
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Minimize  /  =  tf 

Subject  to  v,  =  y4 

y2  =  y5 
y3  =  y6 

y4  =  U\ 

%  =  u2 
y6  =  ui 
y2  =  Jio 
.Vs  =  i 

%  =  yn  .  (4.16) 

-Vio  =  -2^7  +  y8  -  2-VlO  +  Vu  +  “l 

ji  1  =  y2  -  2  +  .To  -  2yt  i  +  yl2  +  U2 

yl2  =  y&-y9+yn-yn+u3 
(yi{t0),y3{t0),y5{t0),y6{t0))  =  {-io,5, o,o.i) 

(.V2^/)^3^/)^5(^)^6('1/))  =  (15  ,-15,o.i,o) 

t0  =  o 

|yj<  0.75,/  =  4,5,6 
\u\  <  1.5,/  =  1,2,3 

The  dynamics  in  (4.16)  are  such  that  constraining  the  total  motion  of  the 
system  is  simple.  However,  they  do  not  readily  provide  modal  characteristics. 
The  dynamics  are  expressed  in  physical  space,  so  the  modes  cannot  be 
constrained  directly.  The  dynamics  must  be  uncoupled  by  transforming  them  to 
modal  space,  thereby  enabling  effective  mode  constraints. 

B.  UNCOUPLING  THE  3-DOF  SYSTEM 

The  system  as  constructed  exhibits  proportional  damping.  Therefore,  the 
undamped  mode  shape  vectors  can  be  used  to  diagonalize  the  mass,  damping, 
and  stiffness  matrices  to  produce  the  modal  matrices,  and  thus  the  system  can 
easily  be  uncoupled  into  three  single-DOF  systems  of  the  form 

m,.qr  +  crc/r  +  k,q,  =  f'  >  (4-17) 
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where  r  is  the  mode  number.  Recall  from  Chapter  II  that  f' =  Vrf .  The  transfer 
functions  in  the  Laplace  domain  can  then  be  found  easily.  Assuming 

q(0)  =  q(0)  =  0 , 

r  k  F' 

Qrs2  +  ^rQrs  +  ^rQr=^r 
mr  mr  mr 

]_  (4.18) 

Q r(s) = Hr(s)F;(S)  =  —  m'  -  f;c.v) 

s 2  +  ^-s  H — — 

mr  mr 

These  transfer  functions  represent  displacement  in  modal  coordinates  due  to  the 
transformed  forces,  and  they  are  much  easier  to  derive  and  manipulate  than  the 
transfer  function  matrix  in  (3.2).  It  is  important  to  distinguish  the  difference 
between  an  uncoupled  transfer  function  in  modal  (principal)  coordinates  versus  a 
coupled  transfer  function  in  physical  coordinates. 

The  uncoupled  equations  of  motion  from  (4.17)  are 

f 

+0.19814,  +0.1981*7!  =  — — — 

f' 

q2  +  l.5550q2  +  l.5550q2=  Jl  .  (4.19) 

1 .8412 

f' 

q3  + 1 .5550*73  + 1 .5550*73  = - - — 

2.8629 

These  equations  represent  uncoupled  modal  masses  as  in  Figure  22.  The 
optimal  control  problem  definition  can  now  be  redefined  in  terms  of  modal 
coordinates.  Let  the  state  vector  be  defined  as 


y  = 


X\  X2  X3  Xf  X2  X3  <7,  q2  q3  <74  q5  qk 


where 


94  9s  9  6 


9i  92  93 


T 

,  and 


-I  T 


X{  X2  X3  X,  X2  X3  <7i  92  93  94  9s  96 


(4.20) 


(4.21) 
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;  m 

q3  =  ° 


Figure  22.  Uncoupled  3-DOF  mass-spring-damper  system. 

-|  T 

The  control  vector  remains  u=  0  f2  f3  .  However,  the  dynamics  equations 

iT 

for  the  flexible  component  defined  by  %  •••  yl2  must  account  for  the 

transformation  from  physical  coordinates  to  modal  coordinates.  For  example, 
from  (4.19)  we  can  see  that 

^1  =  9^9_0-19814l_0-1981?1  '  (4'22) 

Because  f'=Vnf  +V2lf2  +V3lf3  =  m1(vnal  +V2la2  +  V31a3)  and  cjt  =  q4 ,  the  general 
form  is 

?4  =  ?l(^iA  +  ^2i«2  +  ^3A)-^l?4--^?i  ■  (4-23) 

m,  ml  m, 

Thus,  the  general  form  of  the  optimal  control  problem  in  modal  coordinates  can 
be  defined  as  in  (4.24).  Note  that  there  are  no  constraints  or  terminal  boundary 
conditions  placed  on  the  flexible  motion  in  the  problem  statement.  The  effects  of 
constraining  the  flexibilities  are  discussed  in  a  later  section. 
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Minimize 
Subject  to 


J  =  tf 
y,  =  y  4 
%  =  y3 
y3  =  y6 

u. 

y*  =  — 

ml 

«, 

=  — 

m2 


m3 

y  i  =  yw 

%  =  -V„ 


.(4.24) 


y  9  =  y  12 

1  C  k 

yW=—  (VU«l+y2.«2+V31M3)-^>’lO-^>'7 

m]  mx  m] 

1  c  k 

T.  =  —  (y.2«l  +V22U2  +y32«3)-  -=~yn  -  -^y* 
m2  m2  m2 

1  c  k 

*12  =  —  (y!3Ml  +  V23U2  +  ^3“3)-  ~=~yi2  - 

m3  m3  m3 

(};i(?o).>'2(?o)>>,3(?o).3;4(?o)»>'5(?o)>>'6(?o))  =  (0  -10,5,0,0,0.1) 

(>i  )  .ys  (^ )) = (o,  is  -15,0,0.1,0) 

(O’>io  (O’TlCO’TjlO)^0’0’0’0’0’0) 

|v.|  <0.75,  /  =  4,5,6 
|n,.|  <1.5,  i  =  1,2,3 


C.  SLEW  OPTIMIZATION  WITH  NO  CONSTRAINTS  ON  VIBRATIONS 

The  MATLAB  package  DIDO  [6]  was  used  for  optimizing  and  validating 
the  minimum-time  maneuver.  The  minimum-time  maneuver  is  presented  in 
Figure  23.  For  both  axes,  the  maximum  acceleration  (instead  of  some  fraction  of 
it)  is  supplied  initially  and  at  the  end  of  the  maneuver  for  a  much  shorter  time 
than  in  the  computed  torque  case.  It  is  worth  noting  that  y3  does  not  follow  the 
exact  same  trajectory  as  in  the  computed  torque  slew  in  Figure  18.  The 
acceleration  is  not  supplied  in  a  “bang-off-bang”  method  but  rather  oscillates 
slightly.  Because  y2(t/)-y2(0)>y3(t/)-y3(0),  the  azimuth  axis  becomes  the 

limiting  axis.  As  a  result,  the  elevation  axis  has  additional  freedom  in  the  path  it 

71 


takes  to  the  desired  endpoint.  This  is  why  the  acceleration  applied  to  that  axis  is 
not  uniform.  As  an  analogy,  picture  two  baseballs  being  thrown  simultaneously 
from  the  same  initial  position  with  the  same  constraint  on  maximum  velocity,  but 
with  the  first  baseball’s  target  10  feet  further  than  the  second’s.  The  shortest 
travel  time  for  the  first  baseball  dictates  the  shortest  travel  time  total  for  both 
baseballs  together.  The  second  baseball’s  velocity  does  not  need  to  match  the 
first  baseball’s  velocity  to  reach  its  target  at  the  same  time.  For  the  optimal 
control  problem  in  (4.24),  that  time  is  33.8  seconds,  a  savings  of  4.3  seconds 
(12.7%)  over  the  computed  torque  maneuver. 


Figure  23.  Rigid-body  results  of  optimal  control  problem  defined  in  (4.24) 
(tf  =  33.8) :  (a)  position;  (b)  rate;  (c)  acceleration. 


Figure  24  shows  the  resulting  flexible  motion  in  both  modal  coordinates  (in 
terms  of  q)  and  physical  coordinates  (in  terms  of  x).  Recall  that  x  =  Vq  and  that 
qr  ^  .  Rather,  ql  represents  the  motion  of  all  masses  with  respect  to  mode  1 , 
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while  Xj  is  the  physical  displacement  of  mass  1  only,  and  so  on.  Figure  24  also 
illustrates  the  relationship  x  =  Vq.  For  example, 

q(5)  =  [  -0.1165  -0.0231  0.0051  J ,  and  the  second  row  of  the  mode  shape 

matrix  is  [  1.8019  0.4450  -1.2470  ].  Thus,  x2 (5)  =  -0.2266,  as  the  plot  shows. 
It  is  also  evident  that  mode  1  dominates  the  total  motion  of  the  system. 


Figure  24.  Flexible-body  results  of  optimal  control  problem  (4.24): 
(a)  modal  coordinates;  (b)  physical  coordinates. 


The  added  freedom  in  the  path  of  y3  (discussed  previously)  allows  control 
designers  to  apply  cost  reduction  measures  to  the  system.  For  example,  it  is 
possible  to  find  a  minimum-time  maneuver  that  also  minimizes  the  energy 
expended  to  control  the  non-limiting  axis.  In  this  case,  the  problem  definition  in 
(4.24)  can  be  altered  to  minimize  the  cost  function 
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(4.25) 


v 

J  =  200 J  u2  dt  +  tf 

0 

subject  to  the  same  dynamics  and  boundary  conditions/constraints.  The  integral 
term  in  the  cost  function  represents  a  running  cost,  in  that  it  is  calculated  and 
summed  at  each  time  step,  and  it  is  weighted  200  times  against  the  endpoint 
cost  of  final  time.  By  doing  this,  the  control  (energy)  is  reduced  while  also 
minimizing  time.  The  new  angular  rates  and  accelerations  for  the  revised  optimal 
slew  are  shown  in  Figure  25.  The  displacement  of  y3  now  more  closely  matches 
the  computed  torque  slew  while  completing  the  maneuver  in  the  same  time  as 
Figure  23.  Due  to  (4.25),  the  maneuver  realizes  additional  savings  in  energy 
expended. 


(b) 


Figure  25.  Rigid-body  results  of  minimum-time  slew  with  running  cost  to 
minimize  control  effort  [tf  =  33.8) :  (a)  position;  (b)  rate; 

(c)  acceleration. 
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The  new  acceleration  profile  ostensibly  affects  vibrations,  shown  in  Figure  26  and 
Figure  27.  In  Figure  26,  mode  1  again  clearly  dominates  the  total  motion  of  the 
system  but  in  a  more  pronounced  way  as  the  spacecraft  minimizes  its  control 
input.  Modes  2  and  3  show  less  oscillation  than  before  (compare  Figure  24  to 
Figure  26). 


Figure  26.  Flexible-body  results  of  optimal  control  problem  with  running 

cost  in  modal  coordinates. 

Figure  27  shows  a  comparison  of  each  mass’s  displacement  (magnitudes  only) 
for  each  optimal  control  solution  (minimum-time  and  minimum-time  plus  minimum 
control  effort).  The  addition  of  a  running  cost  generally  reduces  the  magnitude  of 
displacement  across  all  three  masses,  which  follows  as  the  applied  control  is 
minimized.  Applying  a  running  cost  to  minimize  applied  control  can  be  a  feasible 
method  of  reducing  vibrations  through  a  maneuver  and  is  one  of  the  tools 
available  to  designers  using  the  optimal  control  approach.  Note  that  this  option  is 
not  possible  using  conventional  control  concepts. 
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The  plots  of  the  modal  coordinate  q  illustrate  an  interesting  phenomenon. 
Recall  from  (3.6)  that  the  damped  frequencies  for  each  mode  are 

"  0.4339  0.9749  0.7818  ]  rad/s,  or  [  0.0691  0.1552  0.1244  ]  Hz.  The 

natural  tendency  for  any  system  is  to  oscillate  at  its  natural  frequency  [7],  The 
same  is  true  for  this  system.  In  Figure  26,  for  example,  the  period  of  q2 
decreases  over  time  before  the  last  acceleration  pulse  is  applied,  eventually 
reaching  about  6.5  seconds,  or  roughly  0.154  Hz. 


Xi  X2 


. Minimum  time  only 

- With  running  cost 


Figure  27.  Comparison  of  flexible-body  results  in  physical  space  for 
endpoint  cost  vs.  endpoint  plus  running  costs  (magnitudes  only). 


The  implementation  of  the  spacecraft’s  control  system  also  becomes 
apparent.  The  displacement  of  stays  at  or  very  near  to  zero  (on  a  scale  of 

10  10 ,  see  Figure  23)  throughout  the  rigid-body  simulation  but  its  flexible-body 
counterpart  does  not.  This  vibrational  motion  is  analogous  to  the  control  system 
counteracting  the  effects  of  the  motion  of  the  two  axes  to  keep  the  spacecraft 
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stable  in  its  position.  As  we  constrain  the  control  input,  and  thus  reduce  the 
applied  force  on  y2  and  y3 ,  it  follows  that  the  relative  magnitude  of  ql  decreases 
since  mode  1  is  closely  mapped  to  the  rigid-body  motion  of  the  system  (as 
previously  discussed).  By  reducing  the  force  used  to  slew  masses  2  and  3,  the 
spacecraft’s  control  system  does  not  need  to  work  as  hard  to  counteract  the 
displacements  of  the  other  axes. 

Due  to  the  excitation  of  the  flexible  modes,  the  spacecraft  and  its 
appendages  will  continue  to  vibrate  after  the  slew  has  completed  (see  Figure  25). 
Figure  25  shows  the  behavior  of  the  vibrations  after  the  final  conditions  are  met, 
when  t  >  tf  and  the  payload  (antenna)  must  be  operated.  In  Figure  28,  the 

physical  properties  of  each  mode  take  over  as  the  slew  control  system 
deactivates.  Modes  2  and  3  dampen  quickly  while  mode  1  takes  longer  to  settle. 


Figure  28.  Flexible  spillover  for  time-optimal  slew  maneuver:  (a)  modal 

coordinates;  (b)  physical  coordinates. 


Further,  referring  to  Figure  28,  as  mode  1  settles,  its  frequency  approaches  its 
damped  natural  frequency.  The  bodies  in  each  axis  continue  to  physically 


77 


oscillate  about  equilibrium  but  at  the  same  frequency  (dominated  by  mode  1). 
From  this  result,  it  is  apparent  that  performing  the  slew  maneuver  excites  the 
base  body  in  a  way  that  perturbs  the  motion  of  the  entire  system.  If  the 
perturbation  is  large  enough,  it  may  be  the  case  that  the  operation  following  the 
slew,  (e.g.,  communicating  over  an  antenna  link),  could  be  disrupted.  Therefore, 
it  is  desirable  to  try  to  limit  this  flexible  motion.  This  will  be  the  subject  of  the 
sections  that  follow. 

D.  SETTING  TERMINAL  CONDITIONS  ON  THE  FLEXIBLE  MOTION 

Suppose  it  is  desired  to  extinguish  all  vibrations  at  the  end  of  the 
maneuver.  The  appropriate  terminal  flex  constraints  added  to  the  problem 
statement  in  (4.24)  are 

(<5r](r/),?2(r/),?3(f/),^4(f/),^(f/),?6(r/))  =  (0,0,0,0,0,0)  .  (4.26) 

One  would  expect  that  more  control  or  more  time  is  required  to  cease  flexible 
motion  and  meet  the  desired  end  conditions  for  the  maneuver.  The  minimum¬ 
time  solution  is  shown  in  Figure  29  and  confirms  this  assertion.  The  maneuver 
now  takes  35.14  seconds,  an  increase  of  1.34  seconds  (4%),  as  more  time  is 
required  to  null  the  vibrations.  Similar  behavior  occurs  in  the  elevation  axis  as 
discussed  previously  since  the  spacecraft  must  travel  a  greater  distance  in  the 
azimuthal  direction.  Note  that  in  Figure  29  that  the  trajectories  are  plotted  at  the 
solution  nodes.  They  would  be  smoothed  by  the  correct  interpolation  prior  to 
implementation  on  a  real  system. 


78 


Figure  29.  Optimized  maneuver  with  zero-flex  conditions  at  endstate 
[tf  =  35.14) :  (a)  position;  (b)  rate;  (c)  acceleration. 


Figure  30  plots  the  maneuver  in  modal  and  physical  coordinates.  Note  the 
impulse  on  all  three  axes  at  ~23  seconds;  this  coincides  with  the  acceleration 
impulse  along  the  elevation  axis  in  Figure  29(c)  and  provides  an  illustrative 
example  of  how  the  three  axes  are  coupled.  Force  is  applied  on  only  one  mass, 
yet  this  induces  motion  on  the  other  masses.  In  Figure  30(a),  mode  1  again 
dominates  the  total  motion. 
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Figure  30.  Flexible  spillover  with  terminal  flex  constraint:  (a)  modal 

space;  (b)  physical  space. 


A  running  cost  is  now  added  to  minimize  control.  Applying  the  cost 
function  in  (4.25),  the  optimized  maneuver  (Figure  31)  takes  35.26  seconds,  only 
marginally  slower  than  in  the  previous  case  (0.12  seconds,  0.36%).  This  shows 
that  it  is  possible  to  reshape  the  trajectory  of  the  non-dominant  slew  axis  in  order 
to  reduce  the  control  effort.  Note  that  the  acceleration  profile  for  the  azimuth  axis 
(y2)  is  essentially  the  same  as  before.  However,  the  reduction  of  the  control 
effort  on  the  elevation  axis  has  an  advantageous  impact  on  the  flexible  effects. 
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Time  (sec) 


Figure  31 .  Maneuver  optimized  for  time  and  control  effort  [tf  =  35.26 ) : 

(a)  position;  (b)  rate;  (c)  acceleration. 


Figure  32  illustrates  the  flexible  components  of  the  minimum  time  plus 
minimum  effort  slew.  The  total  magnitude  of  flexible-body  displacement  has  been 
significantly  reduced.  Variable  x3  alone  experiences  roughly  a  50%  decrease  in 
peak  displacement.  With  the  negligible  increase  in  maneuver  time  and  decrease 
in  control  and  vibrations,  this  maneuver  proves  to  be  more  desirable  from  the 
point  of  view  of  system  operation  since  the  load  on  the  base  body  controller  is 
reduced. 
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Figure  32.  Flexible  effects  from  the  maneuver  in  Figure  31 :  (a)  modal 

space;  (b)  physical  space. 


An  interesting  phenomenon  occurs  in  these  two  maneuvers  where  the 
system  uses  the  additional  freedom  in  the  motion  of  the  elevation  axis  to  meet 
the  desired  endstate  while  maintaining  steady  acceleration  on  the  azimuth  axis. 
The  controller  excites  the  modes  of  the  system,  by  applying  force  only  on  m3,  in 
such  a  way  that  its  modal  properties  work  to  reduce  the  effects  of  flexibility  in  the 
system.  This  non-intuitive  aspect  is  only  possible  by  using  optimal  control  to 
design  the  maneuver  trajectories. 

To  further  illustrate  the  point,  assume  that  the  initial  conditions  are 
changed  to 

(^o)’3^5 (?o)^6(^o))  =  (°  -10,10,0,0,0.1)  (4.27) 

so  that  both  axes  now  must  move  the  same  distance,  such  that 
.^2  (r/) —  (°) =  -v3  (r/) —  -Ts (°)  ’  thus  removing  freedom  previously  present  in  the 

elevation  axis.  The  results  are  telling;  the  minimum-time  maneuver  (no  running 
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cost)  is  0.89  seconds  (2.5%)  slower  (see  Figure  33),  and  the  controller  must  now 
torque  both  axes  before  the  end  of  the  maneuver  in  order  to  suppress  the 
flexibilities.  A  similar  effect  is  observed  for  the  minimum  time  and  effort  problem 
(see  Figure  35).  The  vibration  profiles  for  each  maneuver  (see  Figure  34  and 
Figure  36)  are  similar,  further  emphasizing  the  role  a  non-dominant  axis  can  play 
in  accommodating  flexible  effects. 


(b) 


Figure  33.  Minimum-time  maneuver  for  y2(/L/)-y2(0)  =  y3(r/)-y3(0) 
(tf  =  36.03) :  (a)  position;  (b)  rate;  (c)  acceleration. 
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Flexible  displacement 


(b) 

Time  (sec) 


Figure  34.  Flexible  displacement  from  minimum-time  maneuver  of 
Figure  33:  (a)  modal  space;  (b)  physical  space. 
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(b) 


Figure  35.  Maneuver  optimized  for  time  and  control  effort  with 

y2(^/)-y2(0)  =  ^3(^/)-y3(0)  =  36.13) :  (a)  position;  (b)  rate; 

(c)  acceleration. 
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Figure  36.  Flexible  displacement  from  maneuver  of  Figure  35:  (a)  modal 

space;  (b)  physical  space. 


E.  BOUNDS  ON  FLEXIBILITY  DURING  SLEW 

A  less  stringent  case  than  designing  for  zero  spillover  at  the  end  of  a 
maneuver  is  to  set  a  bound  on  the  flexible  displacement  during  the  maneuver.  In 
other  words,  maintain  flexible  motion  within  a  desired  range.  Perhaps  a  designer 
is  willing  to  trade  flexible  effects  for  slew  time  or  control  effort.  One  cannot  simply 
set  nonzero  bounds  on  q(r)  and  expect  them  to  translate  into  the  physical  world. 

As  discussed  previously,  the  transformation  x  =  Vq  must  be  applied.  For 
example,  suppose  it  is  desired  to  bound  vibrations  in  physical  space  for  all  axes 
to  the  interval  [-0.3,  0.3]  degrees  with  no  constraint  on  x.  The  transformation  for 
the  upper  and  lower  bounds  yield,  for  the  example  problem  of  this  chapter, 

qupper  =  [  0.1629  0.1048  0.0323  J  (4.28) 
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(4.29) 


qhwer  =  [  -0.1629  -0.1048  -0.0323  J  . 

However,  the  simple  box  bounds  given  in  (4.28)  and  (4.29)  do  not  properly 
constrain  q,  resulting  in  potential  violations  of  the  constraint  on  x  (see 
Figure  37). 


Time  (sec) 

Figure  37.  Flexible  displacement  for  a  slew  with  constraints  on  x  applied 
via  x  =  Vq :  (a)  modal  constraints  satisfied;  (b)  physical 
constraints  violated. 

The  constraint  violation  in  Figure  37  occurs  because  each  modal  coordinate 
depends  on  the  values  of  all  three  physical  coordinates  (and  vice-versa  when 
transforming  the  other  direction).  For  a  simple  explanation,  consider  a  2-DOF 
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system  with  eigenvector  matrix  V  = 


.  A  plot  of  the  box  bounds  on  x  is 


a  b 
c  d 

shown  in  Figure  38  along  with  corresponding  bounds  on  q  =  V  'x . 


“Naive”  bounds 


Figure  38.  Example  constraints  on  x  transformed  to  q  via  q  =  V  'x. 


The  area  in  the  figure  outlined  in  solid  black  represents  all  of  the  values  of  q1  and 
q2  in  modal  space  that  satisfy  the  constraints  in  physical  space.  However,  the 
transformation  of  the  bounds  on  x  directly  into  modal  space  produces  the  large 
rectangle  in  dashed  red  (naive  bounds).  In  this  case,  there  are  points  in  modal 
space  that  satisfy  the  constraints  on  x  but  not  on  q.  Plots  of  q  transformed  to  x 
would  show  the  same  behavior.  As  seen  in  Figure  37,  a  transformation  of  bounds 
in  physical  space  to  modal  space  and  back  results  in  violation  of  the  constraints. 

One  approach  to  overcome  this  is  illustrated  by  the  small  rectangle  in 
dashed  red  in  Figure  38.  Instead  of  constraining  x1  and  x2  simultaneously, 
perform  a  transformation  on  each  individually  and  use  the  minimum  absolute 
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values  of  q  to  form  a  conservative  constraint  boundary  as  in  Figure  38.  For  the 
example  3-DOF  system,  the  modal  constraints  [qx^  for  each  axis  are 

qXi  =  V“'  [  0.3  0  0  J=[  0.0323  0.1629  0.1048  ]' 

qX2  =  [  0  0.3  OH  0.0582  0.0725  -0.1307  J  ■  (4.30) 

qx_  =  V“'  [  0  0  0.3  J  =[  0.0725  -0.1307  0.0582  ]' 

The  resulting  constraints  on  each  modal  coordinate  are 

q„,-,«=±[  0.0323  0.0725  0.0582  ]'  (4.31) 

This  approach  creates  conservative  bounds  in  modal  space  but  ensures  that  the 
maneuver  stays  within  desired  bounds  in  the  physical  space.  Interestingly,  when 
optimizing  the  maneuver  for  both  time  and  control  effort,  the  maneuver  looks 
similar  to  the  ones  generated  using  bang-off-bang  control  (see  Figure  39).  The 
maneuver  time  is  34.50  seconds.  The  induced  vibrations  (see  Figure  40)  are 
observed  to  ride  the  constraints  indicating  the  conservatism  inherent  to  the 
maneuver. 
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(C) 

Time  (sec) 

Figure  39.  Maneuver  optimized  for  time  and  control  effort  with 
conservative  constraints  on  q  [tf  =  34.5o) :  (a)  position;  (b)  rate; 

(c)  acceleration. 
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Figure  40.  Flexible  displacement  for  the  maneuver  of  Figure  39: 
(a)  modal  space;  (b)  physical  space. 


As  an  alternative,  a  path  constraint  can  be  used  to  implement  equations 
that  describe  the  true  modal  constraints  shown  in  Figure  38.  In  DIDO,  the 
constraints  are  applied  in  physical  space  by  transformation  of  the  current  modal 
state,  as 


p  =  V 


yi  y%  y9 


(4.32) 


Using  (4.32),  the  constraints  in  physical  space  are  applied  to  the  transformation 
of  the  modal  coordinates  at  each  node.  The  optimal  maneuver  with  the  path 
constraints  is  shown  in  Figure  41  and  Figure  42.  Note  that  the  execution  time  is 
close  to  that  of  the  previous  case,  and  additional  torque  must  be  applied  to  each 
axis  at  certain  instants  in  time  to  avoid  violating  the  constraint  placed  on  the 
amplitude  of  the  vibrations  (compare  Figure  39(c)  with  Figure  41(c)).  Recall  that 
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the  bound  in  physical  space  vibrations  is  -0.3  <x  <0.3;  the  path  constraint 
properly  enforced  this  bound  as  shown  in  Figure  42(b). 


(b) 


Figure  41 .  Maneuver  optimized  for  time  and  control  effort,  using  a  path 
constraint  to  suppress  vibrations  (tf  =  34.54) :  (a)  position; 

(b)  rate;  (c)  acceleration. 
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Figure  42.  Flexible  displacement  for  maneuver  in  Figure  41 :  (a)  modal 

space;  (b)  physical  space. 


These  slew  maneuvers  are  compared  to  the  ones  from  Chapter  III  in 
Table  1.  Not  only  is  it  possible  to  constrain  flexible  motion,  but  slew  times 
decrease  when  using  optimal  control  rather  than  increase  when  using  computed 
torque  control.  The  table  highlights  the  fact  that  optimal  control  provides  the 
designer  with  multiple  options,  each  reducing  slew  times  compared  to  traditional 
control  methods. 
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Table  1 .  Comparison  of  slew  maneuvers  using  computed  torque 

control  and  optimal  control. 


Computed  torque 
control 

Optimal  control 

Time 

savings 

h 

(sec) 

Max  amplitude 
(deg) 

(sec) 

Max  amplitude 
(deg) 

Relative  to 
computed 
torque 

During 

slew 

Spillover 

During 

slew 

Spillover 

Unconstrained  flexibility 

38.1 

0.143 

0.118 

33.8 

0.392 

0.266 

1 1 .3% 

Zero  terminal  flex 

Cannot  execute 

35.14 

1.04 

0 

7.7% 

Zero 

terminal 

flex 

Conservative 

bounds 

Cannot  execute 

34.50 

0.15 

0 

9.5% 

Path 

constraint 

Cannot  execute 

34.54 

0.30 

0 

9.3% 

F.  SUMMARY 

In  this  chapter,  the  antenna  slew  problem  was  reformulated  in  order  to 
apply  optimal  control  theory  to  determine  the  control  effort  required  to  minimize 
slew  time  while  simultaneously  managing  structural  flexibility.  It  was  shown  that  it 
is  possible  to  apply  torque  in  such  a  manner  that  the  flexible  properties  of  the 
system  could  be  used  advantageously  to  control  vibrational  effects.  It  was  also 
shown  that  optimal  control  could  be  utilized  to  minimize  the  control  effort  without 
significantly  increasing  slew  time  or  flexible  motion.  This  effect  is  particularly 
useful  in  instances  where  one  link  has  a  greater  distance  to  travel  than  the  other. 

This  chapter  additionally  illustrated  the  value  of  working  in  modal  space 
versus  physical  space.  Modal  space  readily  produces  the  natural  properties  of 
the  system  by  mode  and  allows  control  system  designers  to  understand  the 
responses  of  a  flexible  system  to  force  inputs.  The  equations  of  motion  in  modal 
space  represent  single-DOF  systems  that  uncouple  the  physical  space 
equations,  making  analysis  and  control  of  the  system  dynamics  much  easier. 

The  challenges  of  defining  the  constraints  on  vibrations  in  physical  and 
modal  spaces  were  discussed.  The  key  difficulty  is  preserving  the  relationship 
x  =  Vq.  Conservative  bounds  in  modal  space  ensure  that  physical  constraints 
are  met,  but  the  system  must  exert  more  effort  than  is  required  to  stay  within 
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those  bounds.  Naive  bounds  on  system  behavior  in  modal  space  may  cause 
violation  of  the  bounds  in  physical  space.  Therefore,  a  path  constraint  that  uses 
the  relationship  x  =  Vq  must  be  used  to  ensure  the  correct  physical  constraints 
are  met. 

The  optimal  control  solutions  presented  in  this  chapter  greatly  reduced 
flexible  effects  while  simultaneously  decreasing  slew  time.  While  computed 
torque  control  from  Chapter  III  required  constant  tweaking  of  the  system  to 
achieve  desired  results,  the  optimal  control  method  allows  for  easy  constraint 
definition  as  long  as  the  problem  is  posed  correctly.  The  optimal  control  approach 
will  be  applied  to  a  more  realistic  antenna  system  in  the  next  chapter. 
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V.  OPTIMAL  CONTROL  OF  A  DOUBLE- 
GIMBAL  ANTENNA  WITH  FLEXIBLE  JOINTS 


This  chapter  introduces  a  double-gimbal  mechanism  (DGM)  to  add  to  the 
system  from  Chapter  IV  to  provide  a  more  realistic  model  of  an  actual  gimbaled 
spacecraft  antenna.  The  equations  of  motion  of  this  system  are  significantly  more 
complex  than  those  of  the  simple  models  in  previous  chapters.  The  rotational 
dynamics  introduce  additional  nonlinear  terms  that  are  dependent  on  the  current 
state  of  the  system.  Building  on  the  developments  of  the  previous  chapter,  this 
chapter  compares  the  results  from  computed  torque  control  and  optimal  control 
of  the  nonlinear  system.  The  goal  of  this  chapter  is  to  show  that  modal  analysis 
and  optimal  control  can  be  applied  effectively  to  nonlinear  systems. 

A.  UPDATING  THE  MODEL 

Figure  43  shows  a  DGM  that  is  now  applied  to  the  model.  The  antenna  is 
attached  to  the  spacecraft  body  through  this  two-axis  gimbal. 


Figure  43.  Diagram  of  two-axis  gimbal  (body  1  represents  azimuth,  body 

2  represents  elevation),  from  [5], 
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The  equations  of  motion  for  this  rigid-body  system  are  derived  in  [5]  as 


/i^  +  /2xrsin202+/2::cos202 

0 

0. 

0 

1 

Si 

CN 

1 

CN 

_ i 

1 

Si 

1 _ 

2(/2«-/2^)^Asin02cos02 
(Aa-AxcK8  m02cos02 


(5.1) 


In  (5.1),  Oy  and  02  represent  the  angular  displacement  of  the  azimuth  and 
elevation,  respectively.  This  model  does  not  account  for  the  motion  of  the 
spacecraft  body  and  ignores  flexibility.  For  example,  the  motion  of  the  azimuth 
joint  coincides  exactly  with  the  appendage  connected  to  that  joint. 

The  model  of  (5.1)  is  now  extended  to  include  flexible  effects.  Figure  44 
shows  a  schematic  of  a  flexible  joint.  The  spring  with  stiffness  constant  kt 
represents  joint  flexibility  and  allows  elastic  displacement  between  the  rotor 
(mounted  on  the  spacecraft)  and  the  rotary  link.  The  dashpot  with  damping 
constant  c;  represents  the  torsional  damping  due  to  rotational  friction,  structural 
properties,  and  external  mechanisms  (in  other  words,  anything  that  dissipates 
energy  during  motion  instead  of  storing  it). 


Figure  44.  Representation  of  an  elastic  joint,  after  [14]. 
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In  Figure  44,  y2i  is  the  angular  displacement  of  axis  i,  and  y2M  is  the  angular 
displacement  of  the  actuator.  These  displacements  are  related  to  the  link  and 
rotor  angles  as 


y2i  =  ei,i  =  l,...,n 


Yu- 


where  n  is  the  number  of  axes  in  the  system  (in  this  case,  n  =  2)  and  G(  is  the 
gear  ratio  of  the  actuator.  For  the  purposes  of  this  model,  it  is  assumed  that  the 
system  consists  of  a  direct  drive  between  the  actuator  and  shaft,  so  there  are  no 
gears  between  the  actuator  and  rotor.  As  a  result,  from  (5.2),  y2i_1  =0,-.  The 

elastic  displacement  in  the  joint,  therefore,  is  y2i  -Yu- 1- 

r  ~\T  r  ~\T 

Let  kt  =  k,  c(.  =  c,  y,=[  y2  y4  J  =[  0,  02  J  ,  and 

-|  T  r  -|  T 

y2  =  7,  y3  =  0i  02  .  The  inertia  tensor  of  the  rotor  is  given  as 

(I,0,0,),=diag[  l,xx  liyy  lizz  ],  and  D(y,)  gives  the  configuration-dependent 

inertia  of  the  rigid-body.  Using  these  definitions,  the  kinetic  energy  for  the  system 
may  be  computed  as 

KE  =  \y\D(yl)y,  +  ^yT2irotoj2  (5.3) 

where  moment  of  inertia  Jrolor  =  diag  Gf (lrotor)la  •••  G;(lrotor)nzz  [14],  The 

potential  energy,  a  function  of  the  angular  displacements  of  the  axes  and  the 
elastic  displacement,  is 

PE  =  P1(yi)  +  P2(y1-y2)  (5.4) 

where  function  P,  is  the  potential  energy  of  the  rigid-body  system.  Function  P2 
depends  on  the  elastic  displacement  and  is  defined  by 
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(5.5) 


1  T 

(yi-y2)  ■ 

The  Lagrangian  may  be  constructed  as 

L  =  KE  —  PE  =  !(yfZ>(y,)yi  +  y2 irotory2  -  2 P1  (y, )  -  k( y,  -  y2 f  (y,  -  y2 ))  (5.6) 

Damping  results  in  energy  dissipation  rather  than  energy  storage.  A  typical 
model  of  this  is  Rayleigh  dissipation  given  as  [15] 

Edis  =  \c(  y,-y2)2  (5-7) 

In  the  Lagrangian  of  (5.6),  dissipation  may  be  added  to  the  total  energy,  and  thus 
the  Euler-Lagrange  equation  becomes 


'  dL]  c>L  |  dEdis  _ 
v  3y,  J  3y;  dy 


The  generic  equations  of  motion  may  now  be  derived  from  (5.8)  as 

7>yi  +  f^y,-^y[-^-y,+^-]  +  c(y1-y2)  +  /:(y1-y2)  =  o  (5.9) 

l  2  dy,  dy, ) 

J,otory2-c(y1-y2)-^(y1-y2)  =  u  .  (5.10) 

Equation  (5.9)  governs  the  motion  of  the  rigid  links,  while  (5.10)  governs  the 
motion  of  the  rotors.  The  second  term  of  (5.9)  in  parentheses  represents  the 
Coriolis,  gravitational,  and  centripetal  forces  acting  on  the  system.  The  new 
r  iT 

control  vector  u=  t,  t2  now  acts  on  (5.10);  its  effects  are  felt  in  (5.9)  via 
transmission  through  the  flexible  joint. 
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The  new  equations  for  the  flexible  DGM  are  thus 


A:  +  4v>n202+4;cos2<92 

0 

0, 
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where  the  link  torques  in  (5.1)  are  replaced  by  the  torque  transmitted  via  terms 

-c(0i-0,)-£(fli-0i)  and  -c(02-<j)2)-k(e2-<p2). 


B.  SIMULATION  OF  THE  COUPLED  MOTION 

Following  from  the  example  from  [5],  let  /j  =  diag[  4  4  6  ], 

/2  =  diag [  224],  (lrotor).  =  diag[  111],  and  c  =  k  =  1  with  appropriate 
units.  The  equations  of  motion  may  thus  be  rewritten  as 


6  +  2  sin2  01  +  4  cos2  02  0 

0'i 

+ 

-40,02  sin(?2  cos  02 
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2 02  sin  d2  cos  6, 

(5.13) 
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A  simple  PD  controller  with  closed  loop  feedback  was  used  to  simulate  a 
conventional  slew  of  the  flexible  system.  The  following  control  law  was  used. 

r  =  Kl,{eJ-e)+KjaJ-a)  (5.15) 
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where  Od  and  cod  are  the  desired  angular  displacements  and  angular  rates  for 
each  axis.  Given  a  natural  frequency  con=  0.6  rad/s  and  damping  ratio  £  =  1  to 
prevent  overshoot,  the  gains  may  be  computed  in  real-time  as  the  inertia 
changes  for  the  system  (this  is  a  version  of  computed  torque  control). 


KP=K 


/uz  +  /2«sin202  +/2^cos202 


l2yy 


Kd  =  2f«)„ 


/i=+/2>.sin' -e,+/2=cos  e2 


l2yy 


(5.16) 


The  slew  was  implemented  with  the  following  boundary  conditions  and 
constraints: 

(#i  (*o  )A  (f0)’®2  (?o ))  =  (40°, 20°, 0,0) 
(6>1(^),6>2(f/),«I(^),«2(^))  =  (0,40°,0,0)  .  (5.17) 

(r0,r/)  =  (0,20)s 

The  damping  and  stiffness  values  affect  how  much  (or  little)  actuation 
(control  torque)  is  needed  to  slew  each  axis  to  the  correct  position.  If  the  stiffness 
value  k ;  — >  °° ,  then  the  elastic  displacement  9i  -  0( .  0 .  The  joint  now  behaves  as 
a  rigid  body.  A  slew  for  the  rigid-body  system  is  shown  in  Figure  45.  Note  the 
similarities  in  the  torque  profiles  for  the  azimuth  axis  for  this  system  and  the  rigid- 
body  systems  studied  earlier  in  Chapter  III.  More  effort  is  required  to  control  the 
azimuth  axis  versus  the  elevation  axis  because  of  the  time-dependent  nonlinear 
inertia  terms. 
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Figure  45.  Computed  torque  control  of  flexible  DGM  with  k  »  0  : 

(a)  angular  displacement;  (b)  angular  rates;  (c)  control  torques. 

Figure  46  shows  the  angular  displacements  and  rates  of  each  axis  and  the 
torques  applied  to  the  rotors  for  k  =  1,  a  much  more  flexible  system.  The  inertia 
term  for  the  azimuth  axis  changes  as  the  elevation  angle  changes,  affecting  the 
gains  and  torques  needed  to  control  the  system.  Note  how  the  controller 
responds  to  the  changing  inertia  in  the  azimuth  axis  (see  Figure  46(c)).  As  noted 
previously,  it  takes  much  more  effort  to  control  the  azimuth  axis  than  the 
elevation  axis,  but  the  torque  profiles  follow  the  same  basic  shapes  as  in 
Figure  45.  Due  to  the  elasticity  in  the  system,  the  controller  must  continue  to 
work  even  after  tf  in  order  to  damp  out  the  vibrational  effects.  The  large  flexible 

displacements  after  the  slew  (post-maneuver  spillover)  could  present  a  problem 
for  antenna  operations. 
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Figure  46.  Computed  torque  control  of  flexible  DGM  with  k  =  1 : 

(a)  angular  displacement;  (b)  angular  rates;  (c)  control  torques. 


Figure  47  illustrates  the  relationship  between  the  axis  angle  9  and  the 
rotor  angle  (j> .  The  elasticity  in  the  joint  causes  the  axis  to  lag  the  rotor  as  the 
actuator  torques  the  rotor.  The  rotor  must  displace  a  certain  amount  before  the 
axis  follows  suit.  As  stiffness  increases,  the  time  histories  of  each  angle 
approach  the  same  values  as  the  lag  decreases. 
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Figure  47.  Time  history  of  rotor,  (/) ,  and  link,  9 ,  motion  for  a  conventional 
slew  with  k  =  1 :  (a)  angular  displacements;  (b)  angular  rates. 

Note  that  in  using  (5.15)  it  is  assumed  that  the  actual  position  of  the 
gimbal  links  can  be  measured.  An  alternate  control  law  could  be  implemented  if 
the  positions  must  be  inferred  at  the  rotor.  In  this  case,  the  torque  law  would  be 

*  =  Kp{Od-<l>)  +  Kd(G>d-i>)  ■  (5-18) 

Due  to  the  elasticity  and  coupling  in  the  system,  the  slew  maneuver  using  (5.18) 
takes  significantly  more  time  to  reach  end  state  (see  Figure  48).  Figure  48 
highlights  one  of  the  challenges  of  feedback  control  of  flexible  systems:  it  makes 
a  difference  where  measurements  are  taken. 
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Figure  48.  Computed  torque  control  of  a  flexible  DGM  using  control  law 
(5.18):  (a)  angular  displacement;  (b)  angular  rates;  (c)  control 

torques. 


C.  OPTIMIZING  THE  MANEUVER 


Because  the  flexible  DGM  system  contains  nonlinear  terms,  the  state- 
space  representation  of  the  equations  of  motion  takes  on  a  different  form. 
Nonetheless,  these  new  dynamics  can  be  used  to  build  an  optimal  control 
problem  to  reduce  the  slew  time.  From  (5.13)  and  (5.14),  the  state  vector  can  be 


assembled  asy  = 


0, 


02  col 


(02 


01  02  01  02 


T 

,  with  rotor  control  vector 


u  = 


Ti 


*2 


T 

.  For  a  minimum-time  slew,  the  resulting  optimal  control  problem 


is 
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Minimize 
Subject  to 


y  i  =  3b 

^2=^4 

•  _  4,v3y4  sin  v2  cosy,  -  y3  +  y7  -  y,  +  y5 
6  +  2  sin2  y,  +  4  cos2  y, 

...  _  ~2y\ sin  y2  cos  y2  -  y4  +  y8  -  y2  +  y6 
2 

j5  =^7 


^6  =  3^8 

%  =  Ttl+y3-y7+y1-y5 
%  =  ^2+y^-ys+y2-y6 
(yMMtv)’yMMto))=(40°’200’0’0) 

(ys  K  )^6  (t0 )  ,yn  {t0 )  ,y8  ( t0 ))  =  (40°,20°,o,o) 

(>i  {tf)^2  (^)»y3(^)>y4  (*/))  =  (o,4o°, o,o) 

(y5  (tf)^6  {tf)’y i  (f/)>y8  (tf))  =  (o,4o°, o,o) 

|.y,|<irad/,/  =  3,4 
\u\  <1.5  N-m,  i  =  1,2 


(5.19) 


The  challenge  in  working  with  these  dynamics,  given  in  physical  space,  is 
that  the  motion  is  highly  coupled,  so  it  is  difficult  to  constrain  the  flexible  effects. 
To  do  this  correctly,  the  dynamics  should  be  given  in  modal  space.  However,  we 
cannot  develop  the  modal  model  of  this  system  in  its  current  form  due  to  the 
presence  of  the  nonlinear  terms.  Through  arithmetic  manipulation,  it  is  possible, 
however,  to  treat  the  nonlinear  terms  as  a  fictitious  torque  input, 


T  = 

4y3y4  siny2  cosy2  -  2zs  sin2  y2  -  4 z5  cos2  y2 

t2 

-2y2  siny,  cosy. 

This  step  uncouples  the  motion  of  0,  and  02  and  provides  a  linear  system  of  the 
familiar  form 


M,0  +  C(e-0)  +  K(0-0)  =  f 

M20-C(0-0)-K(0-0)  =  t 
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(5.21) 


Using  (5.21),  modal  analysis  can  now  be  easily  performed. 

Modal  analysis  of  the  unforced  system  produces  the  following  symmetric 
FRF  matrix: 


H(S)  = 


52  +  5  +  l 

6s4  +  7s3  +  7s2 
0 

5  +  1 

6s4  +  ls3  +  ls2 
0 


0 

52  +  5  +  1 
254  +  353  +  352 

0 

5  +  1 

254  +  353  +  352 


5  +  1 

654  +  753  +  752 
0 

652  +  5  +  1 
654  +  753  +  752 

0 


0 

5  +  1 

254  +  353  +  352 
0 

25  +5  +  1 
254  +  353  +  352 


(5.22) 


Referring  to  (5.22),  the  system  has  four  degrees  of  freedom,  yet  each 
column  of  the  FRF  matrix  includes  only  two  poles.  The  FRF  matrix  shows 
that  torque  applied  to  one  axis  has  no  effect  on  the  other  (recall  that 
X2  =  H2lFx  + H22F2  + H23F3  + H2AFA).  This  makes  sense  since  in  the  modal  model, 
the  motion  of  each  axis  depends  only  on  the  torques  applied  to  that  axis.  Indeed, 
the  coupling  between  the  two  axes  results  from  the  nonlinear  terms,  which  have 
been  relabeled  as  fictitious  torques. 

By  removing  the  nonlinear  terms  it  is  possible  to  analyze  the  motion  of  two 
independent  systems  assuming  the  simple  case  where  f,  =  f2  =  0  (the  next 
section  discusses  how  to  recover  the  full  nonlinear  model).  The  first  system  is 
characterized  by  the  azimuth  terms: 


6 

0 

o, 

+ 

-0i 
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i 

_cd 

i 

_ i 

0 
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0i 

01-0, 

0i-0, 

*1 

The  second  system  is  characterized  by  the  elevation  terms: 


2 
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1 

CN 

:<£>  : 
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(£>. 

to 

1 

to 

_ 1 

+ 

d2  02 
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(5.23) 


(5.24) 
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In  this  context,  an  optimal  control  problem  definition  in  physical  space  may  thus 
be  defined  as  follows: 


Yi  = 

%  = 

...  _  -.^3  +  ^  -  3^1  +  y5 
^5_  6 

.  ->,4+.y»-.y2+>'t 

^4_  2 

y5  =  yi 

%  =  y8 

y7  =  T\+yi-yi+y\-ys  ■  (5.25) 

^8  =^2+^4  -^8+^2  -y6 

(y,  {t0  ),y2  {t0 )  ,y3  {t0 )  >y4  (*o ))  =  (40°,20°,o,o) 

(y5  (*<>  )>y6  {fo  )>y7  {to)>y&  {t0 ))  =  (40°,20°,o,o) 

(^1  (?/  )^2(?/)^3(?/)^4  (^/))  =  (0,40°,0,0) 

(y5  (tf)^6{tf)^i  (tf),y&  (^/))  =  (o,40°, o,o) 

\yt\zl*%,i  =  3A 
\uj\  <1.5  N-m,  i  =  1,2 

Figure  49  shows  the  results  of  the  maneuver  optimization.  However,  because  the 
nonlinear  terms  are  negated,  this  result  is  not  valid  for  the  nonlinear  system. 
Figure  50  illustrates  the  maneuver  that  results  from  implementing  the  control 
history  from  (5.25)  to  control  the  nonlinear  dynamics  of  the  original  problem 
definition  given  by  (5.19).  Clearly,  the  nonlinear  system  does  not  reach  its 
desired  endpoint. 


Minimize 
Subject  to 
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Figure  49.  Optimized  maneuver  for  flexible,  double-axis  gimbal  with 
nonlinear  terms  removed  ( c  =  k  =  1 ):  (a)  link  displacement;  (b) 
link  angular  rates;  (c)  control  torque. 


110 


Figure  50.  Optimal  control  solution  from  Figure  49(c)  applied  to  control 
the  nonlinear  dynamics  of  (5.19):  (a)  link  displacement;  (b)  link 

angular  rates. 


It  is  now  illustrated  how  to  recover  the  original  nonlinear  dynamics  while 
operating  in  the  linear  modal  space.  As  was  shown  previously,  the  nonlinear 
terms  of  the  DGM  may  be  considered  as  a  fictitious  torque  input  by  moving  them 
to  the  right  side  of  the  dynamics  equations.  Modal  analysis  for  this  revised 
system  may  be  carried  out  in  much  the  same  manner  as  in  previous  examples. 
By  treating  the  nonlinear  terms  as  torque  inputs,  the  modal  parameters  are 
effectively  held  constant  throughout  the  maneuver  even  as  the  antenna’s 
configuration  modifies  the  inertia  of  the  system  over  time.  Instead  of  having  to 
update  the  modal  model  at  each  computational  node,  the  dynamics  for  optimal 
control  are  based  on  a  single  nominal  modal  model.  The  mass  matrix 
M  =  diag [  6  2  11  ],  which  produces  a  constant  eigenvector  matrix 
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(5.26) 


10  10 
0  10  1 
10-60 
0  10-2 


Matrix  (5.26)  is  used  to  transform  physical  to  modal  coordinates  (and  vice-versa). 
The  associated  modal  matrices  are 


M  =  diag[  7  3  42  6  | 

C  =  K  =  diag[  0  0  49  9  j 


(5.27) 


Recall  that  for  proportionally-damped  systems  the  eigenvectors  are  the  same  as 
those  of  the  equivalent  undamped  system;  the  modal  matrices  are  scaled  the 
same  as  the  eigenvectors.  These  matrices  yield  the  following  natural  frequencies 
and  critical  damping  ratios: 


co  j 

®3 

co4 


=  co2=  0 
V42 


1.0801 


=  —  =  1 .2247 
2 


Ci  =  C2  =  ° 
r  _  V42  , 
12 
V6 


=  0.5401 

0.6124 


(5.28) 


Note  that  (5.28)  provides  two  rigid  modes  (the  motion  of  the  DGM  links)  and  two 
flexible  modes  representing  the  joint  flex.  The  equations  of  motion  in  modal 
space  are 

1  q  i  =  u[ 

3  O')  W') 

(5.29) 

42g,  +  49g,  +  49^  =  u'3 
6q4  +  9q4  +  9  q4  =  u4 
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where  u'  =  Vru .  The  first  and  second  equations  (for  their  respective  modes) 
result  from  purely  rigid-body  motion. 


A  new  state  vector  y : 


IT 


4l  42  43  44  4l  42  43  44 


is  now  defined 


with  control  u  =  [  f  r  ]  where  f  = 


T,  Z2 


and  r  = 


T  T2 


The 


optimal  control  problem  for  the  nonlinear  system  expressed  in  terms  of  the  fixed 
modal  model  with  fictitious  torque  inputs  is  given  as 


Minimize 
Subject  to 


J  =  tf 

y.  =  >5 
y2  =  y6 
%  =  y7 
y4  =  y8 

y5  =  -(Mi  +  M3) 

y6  =  ^{u2+u4) 

117  7 

y7  =  — u,  — — y7  — y, 
7  42  1  7  3  6 '  6 ' 3 


1  1 


y 8  ^  u2 


y8-^y4 


6  3  2  2 
(Vj  (?0)^2  (?o)^3(?o)^T4(^o))  =  V_1  [ 

(Xs  (?o)^T6  (?o)^7(?o)^V8  (?0))  =  V_1  [ 

(vi(?/),^2(0),^3(0)’T4(?/))=  V“'[ 

(  v5  )  ,  V6  (v  )  ,3^7  (0  )  ^8  (?/  ))  =  V_1 

|  v,|  —  V-1  [  1  1  1  1  Jra%,/  =  5,6,7,8 

M0-t(0I 


40' 

3 

20°  40° 

20"  J 

iT 

0 

0 

0  0  J 

0 

40°  ±180° 

±180° 

0 

0 

+i8o  y 

/sec 

±180/ 

sec 


,(r)-f2(t 

MO  I 


\u4{t 


< 


0 

0 

1.5  N-m 
1.5  N-m 


(5.30) 
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In  solving  (5.30),  the  torque  limits  are  applied  in  modal  coordinates  through  the 
transformation  in  the  dynamics  and  the  angular  rate  limit  at  each  node  as  path 
constraints.  Path  constraints  are  also  used  to  ensure  that  fictitious  torques  f 
admit  the  correct  nonlinear  effect.  It  is  the  presence  of  these  two  path  constraints 
that  allows  the  full  nonlinear  dynamics  to  be  incorporated  with  the  fixed  modal 
model.  This  aspect  represents  a  new  innovation  of  this  thesis. 

As  in  Chapter  IV,  we  first  study  the  behavior  of  the  system  without  any 
constraints  on  flexibility.  Figure  51  shows  the  optimized  maneuver  resulting  from 
this  approach  after  transformation  back  into  physical  coordinates.  The  azimuth 
axis  rotates  the  greater  angular  distance,  requiring  the  greater  control  torque,  but 
the  controller  is  able  to  take  advantage  of  the  added  freedom  in  the  elevation 
axis.  The  torque  profiles  of  Figure  51(c)  are  quite  different  than  those 
corresponding  to  computed  torque  control.  This  allows  the  maneuver  time  to  be 
reduced. 
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Figure  51 .  Optimal  slew  for  modal  system  with  fictitious  torques:  (a)  axes 
displacement;  (b)  axes  angular  rates;  (c)  input  torques. 

Since  we  are  only  concerned  with  the  torques  supplied  by  the  spacecraft 
actuators,  only  z1  and  t2  are  needed  to  implement  the  maneuver.  The  values  of 
the  fictitious  controls  are,  however,  included  for  reference  (see  Figure  52). 
Applying  t,  and  t2  to  control  the  full  nonlinear  dynamic  model  produces  the 
correct  maneuver,  as  shown  in  the  first  part  of  Figure  53  (to  the  left  of  the  vertical 
dashed  line).  Thus,  the  new  optimal  control  problem  formulation  allows  the 
nonlinear  effects  to  be  properly  recovered. 


115 


Figure  52.  Fictitious  control,  f . 


In  solving  the  optimal  control  problem,  the  vibrations  have  not  been 
constrained  in  any  way.  Figure  53  also  shows  the  behavior  of  the  system  after 
the  end  of  the  minimum-time  slew  (to  the  right  of  the  vertical  dashed  line). 
Because  of  the  energy  stored  by  the  “spring”  in  the  joint,  each  rotor  and  axis 
continue  to  spin.  This  maneuver  is  not  feasible  for  operations  since  it  is  required 
for  the  antenna  to  remain  pointed  at  its  target. 
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Figure  53.  Flexible  spillover  for  minimum-time  slew  of  flexible  double 

gimbal:  (a)  axis;  (b)  rotor. 


To  complete  the  maneuver  with  zero  flexible  motion,  the  endpoint 
constraints  in  (5.30)  can  be  modified  to 

(y.  (^ )  (^ )  .y3 (^ )  (^ )) = vr  [  o  40°  o  40°  ] 

(5.31) 

0  0  0  0  ]  =(0,0, 0,0) 

This  ensures  that  the  rotor  and  link  angles  match  at  the  end  of  the  maneuver  so 
that  there  is  no  residual  energy  stored  in  the  spring.  One  would  expect  the  slew 
time  to  increase  because  the  controller  must  now  force  the  total  motion  of  the 
system  to  zero.  Figure  54  illustrates  that  this  conjecture  is  correct,  while  Figure 
55  verifies  that  motion  ceases  at  tf .  Note  also  that  in  the  zero  flex  case  the 

flexible-body  modes  q3  and  q4  are  excited  (see  Figure  56)  because  they  are  not 
otherwise  constrained. 
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Figure  54.  Minimum-time  maneuver  with  endpoint  flex  constraint, 
tf  =  6.496 ,  k  =  1 :  (a)  axes  displacement;  (b)  axes  angular  rates; 

(c)  input  torques. 


(a)  (b) 

Time  (sec)  Time  (sec) 


Figure  55.  Flexible  spillover  for  minimum-time  slew  of  flexible  double 
gimbal  with  terminal  flex  constraints:  (a)  axis;  (b)  rotor. 
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Figure  56.  Displacement  of  system  in  modal  space,  zero  flex  constraints. 


We  can  easily  study  the  effects  of  limiting  flexible-body  motion  now  that 
we  are  operating  in  modal  space.  In  the  context  of  this  problem,  such  a  constraint 
would  limit  the  relative  displacement  between  the  rotor  and  the  link.  To  illustrate 
this,  the  following  constraint  may  be  added: 


1^1  <10° 

kk10° 


(5.32) 


Figure  57  shows  that  due  to  (5.32)  the  maneuver  time  increases  dramatically 
because  the  control  torque  must  be  modulated  to  manipulate  the  relative 
displacement  at  each  joint. 
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Figure  57.  Minimum-time  maneuver  with  flex  constraints  ( k  =  1,  tf  =  1 1 .5 1) 

:  (a)  axes  displacement;  (b)  axes  angular  rates;  (c)  input  torques. 

Figure  58(a)  shows  that  for  this  solution  the  flexible-body  modes  are 
properly  constrained  and  that  the  rigid-body  modes  respond  better  to  the  input 
torques  since  much  less  energy  is  stored  in  the  spring  of  the  joints.  Figure  58(b) 
confirms  that  the  rotor  and  axis  more  closely  track  each  other,  showing  that 
constraining  the  flexible  modes  allows  the  elastic  displacement  to  be  controlled. 
The  final  displacements  of  the  rotors  and  links  are  the  same,  so  the  antenna 
does  not  oscillate  at  the  end  of  the  slew. 
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Figure  58.  Flexible  effects  for  maneuver  with  flex  constraints:  (a)  modal 

space;  (b)  physical  space. 


If  the  endpoint  flex  constraint  on  vibrations  were  removed,  the  antenna 
would  continue  to  move  after  the  controller  shuts  off  (see  Figure  59).  The  motion 
is,  however,  slower  than  that  in  Figure  53  because  of  the  imposed  constraint  on 
the  displacements  of  the  flexible  modes.  Thus,  the  proper  combination  of 
constraints  on  the  flexible  modes  both  during  the  slew  maneuver  and  at  the 
terminal  time  appears  to  be  the  best  way  to  reduce  slew  times  for  flexible 
systems. 
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Figure  59.  Link  angle  for  maneuver  with  constrained  flexibility. 

A  comparison  of  the  slew  maneuvers  for  the  DGM  in  this  chapter  is  shown 
in  Table  2.  The  computed  torque  maneuver  (without  the  nonlinear  terms) 
executes  nearly  five  times  slower  than  the  optimal  control  maneuvers.  The  rigid 
DGM  executes  the  fastest  but  does  not  take  into  account  the  flexible  effects  of 
the  system.  The  optimal  control  solutions  for  the  flexible  DGM  provide  the  most 
realistic  representation  of  the  slew  maneuvers  while  enabling  designers  to  set 
constraints  on  flexible  motion.  Even  with  the  constraints,  these  slew  maneuvers 
execute  faster  than  those  using  conventional  control. 


Table  2.  Comparison  of  slew  times  (in  seconds)  for  a  flexible  DGM. 


Computed 
torque  control 

Optimal 

control 

Time 

Savings 

Unconstrained 

20 

4.5 

77.5% 

Zero  terminal  flex 

Cannot  execute 

6.5 

67.5% 

Zero  terminal  flex 
and  path 
constraint 

Cannot  execute 

11.5 

42.5% 
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D.  SUMMARY 

This  chapter  demonstrated  the  application  of  the  optimal  control  concepts 
developed  in  Chapter  IV  to  a  nonlinear  dynamical  system  having  characteristics 
of  a  real  spacecraft  antenna.  The  effects  of  constraining  modal  displacement 
both  during  the  maneuver  and  at  the  terminal  time  were  discussed.  The  fastest 
slews,  as  expected,  occurred  with  unconstrained  modal  displacement.  These, 
however,  require  control  effort  after  the  slew  to  eliminate  the  spillover  motion  of 
the  antenna.  This  chapter  showed  that  the  most  effective  method  to  constrain 
flexibility  is  to  define  a  path  constraint  and  a  terminal  constraint  on  flexible 
motion.  This  ensures  a  time-optimal  solution  can  be  produced  that  limits 
vibrations  and  ensures  the  antenna  will  remain  pointed  at  the  desired  target  after 
the  slew  is  complete. 

The  key  innovation  in  Chapter  V  and  this  thesis  is  the  separation  of  the 
nonlinear  terms  from  the  model  by  representing  them  as  a  fictitious  torque.  This 
allowed  a  nominal  modal  model  to  be  utilized  for  optimal  control  while  retaining 
the  nonlinear  effects  as  part  of  the  optimal  control  solution.  Using  this  approach, 
the  modal  properties  of  the  system  become  readily  apparent,  and  the  optimal 
control  solution  operates  as  desired  when  applied  to  the  full  nonlinear  dynamical 
system  in  physical  space. 
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VI.  CONCLUSION 


This  thesis  addressed  a  problem  in  satellite  design  and  operations  where 
a  control  system  must  accommodate  the  effects  of  slewing  a  flexible  appendage 
while  maintaining  accurate  orientation  and  minimizing  disturbances  on  the 
spacecraft  bus  and  other  payloads.  The  TDRS,  for  example,  has  limited 
availability  of  its  single  access  antennas  due  to  the  capabilities  of  its  current 
control  systems.  Typical  slew  maneuvers  are  executed  slowly  to  reduce 
disturbances  and  avoid  exciting  flexible  modes,  thereby  reducing  customer 
access  times.  This  thesis  resolved  this  issue  by  determining  optimal  control  effort 
for  rapid  slew  maneuvers  by  taking  advantage  a  satellite  structure’s  modal 
properties  rather  than  discounting  them.  This  thesis  presented  a  new  way  of 
finding  optimal  controls  for  flexible  systems  with  complicated  nonlinear  dynamics. 
Typically,  the  nonlinearities  are  simplified  or  assumed  away  to  make  analysis 
easier.  However,  this  often  yields  inaccurate  results,  especially  for  large  slew 
maneuvers. 

A  simple  linear  mass-spring-damper  system  was  used  to  illustrate  the 
value  of  modal  analysis  to  understanding  the  behavior  of  flexible  structures. 
Modal  analysis  of  a  notional  spacecraft  antenna  was  then  explored  by  using  a 
3-DOF  mass-spring-damper  system  model  where  applied  forces  cause  coupled 
motion  (translation)  of  each  DOF.  Computed  torque  control  was  used  to 
demonstrate  that  flexible  effects  can  only  be  constrained  by  increasing  slew  time, 
which  counters  the  goal  for  control  system  design. 

Pontryagin’s  Principle  was  applied  next  to  the  same  3-DOF  system,  and 
an  optimal  control  problem  was  defined.  Through  modal  analysis  of  the  system, 
the  motion  of  each  DOF  was  uncoupled  in  modal  space.  Flexible  effects  could 
then  be  constrained  in  modal  space,  so  that  the  optimal  control  solution  could  be 
used  to  reduce  slew  times  while  simultaneously  avoiding  excitation  of  the  modes. 
These  results  showed  that  it  is  possible  to  apply  motive  forces  in  a  manner  that 

utilizes  flexible  properties  advantageously. 
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Optimal  control  and  modal  analysis  were  applied  to  a  simplified  model  of 
an  actual  spacecraft  antenna  mounted  on  a  DGM.  The  DGM  was  updated  with 
elastic  joints  to  model  flexibility.  The  nonlinear  terms  in  the  equations  of  motion 
presented  a  problem  for  modal  analysis  in  that  the  modal  model  was  time- 
dependent.  A  novel  approach  was  used  to  recast  the  nonlinearities  as  fictitious 
torque  inputs  to  the  system.  This  approach  removed  the  time-dependency,  yet 
provided  a  nominal  modal  model  that  accurately  described  the  motion  of  the 
original  nonlinear  system.  It  also  enabled  the  application  of  effective  flex 
constraints  for  optimal  control  solutions.  In  the  optimal  control  solution,  the 
controller  used  the  available  torque  to  execute  a  rapid  slew  while  meeting  the 
desired  flex  constraints.  The  maneuver  time  for  the  flexible  DGM  was  decreased 
by  approximately  42%  compared  to  the  conventional  computed  torque  control 
solution.  Thus,  the  throughput  of  an  antenna  system  such  as  the  TDRS  single 
access  antenna  could  theoretically  be  improved  using  the  results  of  this  thesis. 

A.  TRANSITION  TO  PRACTICE 

Transforming  the  dynamics  of  a  system  from  physical  space  to  modal 
space  has  other  advantages.  Realistically,  a  control  system  designer  would  be 
given  the  natural  frequencies  and  damping  ratios  of  a  structure  rather  than 
having  to  derive  them.  If  this  is  the  case,  optimizing  the  control  is  simple.  The 
natural  frequencies  and  damping  ratios  for  each  mode  easily  translate  to  the 
coefficients  of  the  modal  space  dynamics  equations,  as  demonstrated  in 
Chapter  II.  These  coefficients  comprise  the  modal  mass,  damping,  and  stiffness 
matrices  and  are  used  to  determine  the  modal  behavior  of  the  system.  This 
model  can  then  be  incorporated  with  a  nonlinear  dynamics  model  as  was  done  in 
Chapter  V  to  develop  rapid  slew  maneuvers  that  accommodate  flexible  effects 
according  to  the  specification  of  the  mission. 

B.  FUTURE  WORK 

Future  research  can  explore  the  modeling  of  the  behavior  of  systems  with 
hysteretic  and  non-proportional  damping.  Uncoupling  such  systems  is  far  more 
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complicated  than  those  presented  in  this  thesis.  The  simplified  DGM  used  in  this 
thesis  assumed  zero  distance  between  the  azimuth  and  elevation  rotors  and  the 
spacecraft  body;  future  models  can  incorporate  these  offsets  to  study  the  impact 
of  more  complicated  dynamics  on  the  controller’s  ability  to  minimize  vibrations 
and  maintain  accurate  pointing.  This  model  also  assumed  unit  values  for  stiffness 
and  damping  for  both  axes  and  did  not  address  the  impact  of  different  constants 
for  each  axis.  Transition  to  practice  would  require  working  with  the  correct  model 
for  the  flexible  appendage  in  question. 

This  thesis  applied  elasticity  to  the  gimbal  joints  of  a  spacecraft  antenna  to 
model  the  flexibility  of  the  system.  Another  method  would  be  to  model  the 
behavior  of  a  flexible  cantilever  attached  rigidly  to  a  base  structure.  The  results  of 
such  a  model  could  then  be  added  to  those  of  a  model  with  elastic  joints  to 
extend  the  applicability  of  the  results. 
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